#include "macros.h" ! GUGA Data module: ! containing all necessary declarations and basic function definitions for the ! GUGA approach. module guga_data use SystemData, only: nBasis, tSPN, tHPHF, lNoSymmetry, STOT, nEl, & tNoBrillouin, tExactSizeSpace, tUHF, tGUGA, nSpatOrbs_i8, ij_pairs use constants, only: dp, Root2, OverR2, n_int, int_rdm, bn2_, int64 use MemoryManager, only: TagIntType use bit_rep_data, only: BitRep_t, GugaBits use util_mod, only: EnumBase_t, swap, near_zero use indexing_mod, only: fuse_symm_idx, fuse_idx, inv_fuse_symm_idx, inv_fuse_idx use error_handling_neci, only: stop_all implicit none private public :: ExcitationInformation_t, ExcitationType_t, WeightData_t, projE_replica, & getSingleMatrixElement, getDoubleMatrixElement, getMixedFullStop, & orbitalIndex, funA_0_2overR2, minFunA_2_0_overR2, funA_m1_1_overR2, & funa_2_0_overr2, getdoublecontribution, tag_excitations, & tag_tmp_excits, tag_proje_list, funa_3_1_overr2, minfuna_0_2_overr2, & tGUGACore, init_guga_data_procptrs, & excit_type, gen_type, & rdm_ind_bitmask, pos_excit_lvl_bits, pos_excit_type_bits, & n_excit_lvl_bits, n_excit_type_bits, n_excit_info_bits, & RdmContribList_t, n_excit_index_bits, & DistinctDoubleType_t, distinct_doubles, distinct_doubles_arr, & DistinctDouble_t, DistinctDoubleTypeVals_t, & excit_type_arr, valid_excit_type_arr, non_trivial_double_excit_type_arr, & ijkl_Index_t, ijab_Index_t, fused_to_ijkl, fused_to_ijab, fused_to_DistinctDouble ! define types for the probabilistic weights functions used in the ! stochastic excitations generations type :: WeightData_t real(dp) :: F = 0.0_dp real(dp) :: G = 0.0_dp real(dp) :: minus = 0.0_dp real(dp) :: plus = 0.0_dp real(dp) :: zero = 0.0_dp end type WeightData_t type :: GeneratorType_Values_t integer :: & R = 1, & L = -1, & W = 0 end type GeneratorType_Values_t type(GeneratorType_Values_t), parameter :: gen_type = GeneratorType_Values_t() ! if I want to also encode information about the excit_type and the ! excitation level in certain bits of the rdm_ind I think I need to ! create bitmasks to quickly en- and decode this information from the ! integer ! create one general bitmask which indicates the position of the ! excitation type and excitation level information ! and i guess the highest bit is always for the sign.. ! i need 6 bits for the excit-type and 2 for the excitation level ! shift 8 ones to the left (to the most significant bit, except to the ! sign bit (highest bit!) change of mind: lets try using the last bit too! ! and do not use bitmasks for the integers but directly encode them ! and shift them by the correct number of positions in the int_rdm ! and I will use the intrinsic function mvbits() and ibits() to ! perform the en- and decoding ! i need to encode 26 numbers, 2^5 = 32 integer, parameter :: n_excit_type_bits = 5 ! i need 4 values to the excit lvl. 2^2 = 4 integer, parameter :: n_excit_lvl_bits = 2 ! the total number of bits is: integer, parameter :: n_excit_info_bits = n_excit_type_bits + n_excit_lvl_bits ! the position in the bit representation of these information is also ! a helpful quantity: remember bit_size() gives me 64 for a 64bit integer e.g. integer, parameter :: pos_excit_type_bits = int(bit_size(1_int_rdm)) - n_excit_type_bits integer, parameter :: pos_excit_lvl_bits = int(bit_size(1_int_rdm)) - n_excit_info_bits integer :: i_ ! create an integer with as many ones as needed (2**0 + 2**1 + ..) and ! then shift it to the correct position (with the correct amount) integer(int_rdm), parameter :: excitLvl_bitmask = & ishft(int(sum([(2**i_, i_=0, n_excit_lvl_bits - 1)]), int_rdm), & pos_excit_lvl_bits) ! and do the same for the excit_type mask: integer(int_rdm), parameter :: excitType_bitmask = & ishft(int(sum([(2**i_, i_=0, n_excit_type_bits - 1)]), int_rdm), & pos_excit_type_bits) ! and the whole bit-mask is then the or operation of these two integer(int_rdm), parameter :: excitInfo_bitmask = & ior(excitLvl_bitmask, excitType_bitmask) ! and maybe for performance also set a bit-mask of the index part of ! the rdm_ind integer(int_rdm), parameter :: rdm_ind_bitmask = not(excitInfo_bitmask) ! define similar stuff as above for encoding the minimal excit-type ! information into an 64bit integer integer, parameter :: n_excit_index_bits = 8 ! define a type structure to store excitation information between two ! CSFs needed in the matrix element calculation between them ! this may also be used/needed for the excitation generation type, extends(EnumBase_t) :: ExcitationType_t character(30) :: str end type type :: ExcitationTypeValues_t type(ExcitationType_t) :: & invalid = ExcitationType_t(0, 'invalid'), & ! indicate invalid excitation weight = ExcitationType_t(1, 'Weight'), & ! pure weight identifier single = ExcitationType_t(2, 'Single'), & ! all kind of single excitations which dont need much care raising = ExcitationType_t(3, 'Weight + Raising'), & ! weight raising ! W_R(i) ^R(j) ! e_{ii, ij} ! Table 4.3 and Figure 4.4 (0c) lowering = ExcitationType_t(4, 'Weight + Lowering'), & ! weight lowering ! _L(i) W^L(j) ! e_{ji, jj} ! Table 4.3 and Figure 4.4 (0f) non_overlap = ExcitationType_t(5, 'Non Overlap'), & ! non overlap ! e_{ij,kl} or e_{ji,lk} ... ! Table 4.3 and Figure 4.4 (3c0 or 3d0 or 3f0) single_overlap_lowering = ExcitationType_t(6, 'Single Overlap Lowering'), & ! single overlap 2 lowering ! _L(i) ^LL_(j) ^L(k) ! e_{ji,kj} ! Table 4.3 and Figure 4.4 (0e) single_overlap_raising = ExcitationType_t(7, 'Single Overlap Raising'), & ! single overlap 2 raising ! _R(i) ^RR_(j) ^R(k) ! e_{ij,jk} ! Table 4.3 and Figure 4.4 (0b) single_overlap_L_to_R = ExcitationType_t(8, 'Single Overlap L to R'), & ! single overlap lowering into raising ! _L(i) ^LR_(j) ^R(k) ! e_{ji,jk} ! Table 4.3 and Figure 4.4 (1a) single_overlap_R_to_L = ExcitationType_t(9, 'Single Overlap R to L'), & ! single overlap raising into lowering ! _R(i) ^RL_(j) ^L(k) ! e_{ij,kj} ! Table 4.3 and Figure 4.4 (1b) double_lowering = ExcitationType_t(10, 'Double Lowering'), & ! normal double two lowering ! _L(i) -> LL_(j) -> ^LL(a) -> ^L(b) ! e_{ki,lj} / e_{li,kj} ! Table 4.3 and Figure 4.4 (3b) double_raising = ExcitationType_t(11, 'Double Raising'), & ! normal double two raising ! _R(i) -> _RR(j) -> ^RR(k) -> ^R(l) ! e_{jk,il} ! or ! _R(i) -> _RR(j) -> RR^(k) -> ^R(l) ! e_{jl,ik} ! Table 4.3 and Figure 4.4 (3a) double_L_to_R_to_L = ExcitationType_t(12, 'Double L to R to L'), & ! lowering into raising into lowering ! _L(i) -> _RL(a) - > ^RL(j) -> ^L(b) ! e_{jk,li} ! Table 4.3 and Figure 4.4 (3d1) double_R_to_L_to_R = ExcitationType_t(13, 'Double R to L to R'), & ! raising into lowering into raising ! _R(i) -> _LR(a) - > ^LR(j) -> ^R(b) ! Table 4.3 and Figure 4.4 (3c1) double_L_to_R = ExcitationType_t(14, 'Double L to R'), & ! lowering into raising double ! _L(j) -> _RL(a) -> RL^(b) -> ^R(j) ! Table 4.3 and Figure 4.4 (3f1) double_R_to_L = ExcitationType_t(15, 'Double R to L'), & ! raising into lowering double ! _R(j) -> _LR(a) -> LR^(b) -> ^L(j) ! Table 4.3 and Figure 4.4 (3e1) fullstop_lowering = ExcitationType_t(16, 'Fullstop Lowering'), & ! full stop 2 lowering ! _L(i) -> _LL(j) -> ^LL^(a) ! e_{ki,kj} ! Table 4.3 and Figure 4.4 (1d) fullstop_raising = ExcitationType_t(17, 'Fullstop Raising'), & ! full stop 2 raising ! _R(i) -> _RR(j) -> ^RR^(a) ! e_{ik,jk} ! Table 4.3 and Figure 4.4 (1c) fullstop_L_to_R = ExcitationType_t(18, 'Fullstop L to R'), & ! full stop lowering into raising ! e_{jk,ki} ! _L(i) -> _RL(a) -> ^RL^(j) ! Table 4.3 and Figure 4.4 (1e) fullstop_R_to_L = ExcitationType_t(19, 'Fullstop R to L'), & ! full stop raising into lowering ! e_{kj,ik} ! _R(i) -> _LR(a) -> ^LR^(j) ! Table 4.3 and Figure 4.4 (1f) fullstart_lowering = ExcitationType_t(20, 'Fullstart Lowering'), & ! full start 2 lowering ! e_{ji,ki} ! _LL_(i) -> ^LL(a) -> ^L(b) ! Table 4.3 and Figure 4.4 (1h) fullstart_raising = ExcitationType_t(21, 'Fullstart Raising'), & ! full start 2 raising ! e_{ij,ik} ! _RR_(i) -> ^RR(a) -> ^R(b) ! Table 4.3 and Figure 4.4 (1g) fullstart_L_to_R = ExcitationType_t(22, 'Fullstart L to R'), & ! full start lowering into raising ! _RL_(i) -> ^LR(b) -> ^R(j) ! e_{ik,ji} ! Table 4.3 and Figure 4.4 (1j) fullstart_R_to_L = ExcitationType_t(23, 'Fullstart R to L'), & ! full start raising into lowering ! _RL_(i) -> ^RL(j) -> ^L(b) ! e_{ij,ki} ! Table 4.3 and Figure 4.4 (1i) fullstart_stop_alike = ExcitationType_t(24, 'Fullstart to Fullstop Alike'), & ! full start into full stop alike ! _RR_(a) -> ^RR^(i) ! or _LL_(i) > ^LL^(a) ! e_{ij,ij} ! Table 4.3 and Figure 4.4 (2a,2b) fullstart_stop_mixed = ExcitationType_t(25, 'Fullstart to Fullstop Mixed') ! full start into full stop mixed ! _RL_(i) -> _RL_(j) ! e_{ij,ji} ! Table 4.3 and Figure 4.4 (2c) end type interface ExcitationType_t module procedure construct_ExcitationType_t end interface type(ExcitationTypeValues_t), parameter :: excit_type = ExcitationTypeValues_t() type(ExcitationType_t), parameter :: & excit_type_arr(0:25) = & [excit_type%invalid, excit_type%weight, excit_type%single, & excit_type%raising, excit_type%lowering, excit_type%non_overlap, & excit_type%single_overlap_lowering, excit_type%single_overlap_raising, & excit_type%single_overlap_L_to_R, excit_type%single_overlap_R_to_L, & excit_type%double_lowering, excit_type%double_raising, & excit_type%double_L_to_R_to_L, excit_type%double_R_to_L_to_R, & excit_type%double_L_to_R, excit_type%double_R_to_L, & excit_type%fullstop_lowering, excit_type%fullstop_raising, & excit_type%fullstop_L_to_R, excit_type%fullstop_R_to_L, & excit_type%fullstart_lowering, excit_type%fullstart_raising, & excit_type%fullstart_L_to_R, excit_type%fullstart_R_to_L, & excit_type%fullstart_stop_alike, excit_type%fullstart_stop_mixed], & valid_excit_type_arr(1:25) = excit_type_arr(1:25), & non_trivial_double_excit_type_arr(8:25) = excit_type_arr(8:25) type, extends(EnumBase_t) :: DistinctDoubleType_t !! This represents a distinct double excitation (compare Table 4.3 of [@DobrautzThesis2019]) !! A distinct double excitation consists of the two double excitations !! \( e_{ij,kl} \) which are connected by exchange \( e_{il,kj} \). !! !! Because of \( e_{ij,kl} = e_{kl,ij} \) we can order such that \(i \leq k \). !! The matrix element requires \( e_{ij,kl} \) and \( e_{il,kj} \). !! For the notation of distinct double excitation we will order such that \( j \leq l \). !! In the following we assume \( i < j < k < l \). !! !! Example: !! The distinct type ``iijk`` consists of \( e_{ii,jk} \) and \( e_{ik,ji} \). !! The first is a non-overlap double excitation with a single Weight on `i` !! and a single raising excitation from `k` to `j`; !! the latter is a `fullstart_L_to_R` i.e. \(\underline{RL}(i) R\overline{L}(j) \overline{R}(k)\) type(ExcitationType_t) :: repr !! The two subtypes usually do not differ. !! If the two subtypes differ, then it is always a pairing with `non_overlap`. !! For this reason we can choose one representative. !! !! For example: !! !! * `kikj` and `kjki` are both `fullstop_lowering`. !! * `ijkl` is `non_overlap`, the corresponding exchanged excitation !! `ilkj` is `double_R_to_L_to_R`, hence we choose `double_R_to_L_to_R` !! as representation. !! !! @note !! Older code does not clearly distinguish between `DistinctDoubleExcitation_t` !! and `ExcitationType_t`, so the `repr` is particularly useful here. !! @endnote !! type(ExcitationType_t) :: other !! The other Excitation that is not the representative character(2) :: label !! The label in Table 4.3 of [@DobrautzThesis2019] integer :: idx_for_repr(4) = [1, 2, 3, 4] !! Store the permutation for outputting the representation index !! !! If we have for example the `DistinctDoubleType_t` `jilk`, !! then the index of the representation shall refer !! to the `double_L_to_R_to_L` excitation and be `jkli` !! instead of the `non_overlap` part which would be `jilk`. !! The swap is only relevant for distinct doubles with !! a non-overlap excitation. Hence we default-initialize !! to the identity permutation. end type type :: DistinctDoubleTypeVals_t type(DistinctDoubleType_t) :: & invalid = DistinctDoubleType_t(0, excit_type%invalid, excit_type%invalid, ''), & jijk = DistinctDoubleType_t(1, excit_type%single_overlap_L_to_R, excit_type%single_overlap_L_to_R, '1a'), & ijkj = DistinctDoubleType_t(2, excit_type%single_overlap_R_to_L, excit_type%single_overlap_R_to_L, '1b'), & ikjk = DistinctDoubleType_t(3, excit_type%fullstop_raising, excit_type%fullstop_raising, '1c'), & kikj = DistinctDoubleType_t(4, excit_type%fullstop_lowering, excit_type%fullstop_lowering, '1d'), & jikk = DistinctDoubleType_t(5, excit_type%fullstop_L_to_R, excit_type%non_overlap, '1e', [1, 4, 3, 2]), & ijkk = DistinctDoubleType_t(6, excit_type%fullstop_R_to_L, excit_type%non_overlap, '1f', [1, 4, 3, 2]), & ijik = DistinctDoubleType_t(7, excit_type%fullstart_raising, excit_type%fullstart_raising, '1g'), & jiki = DistinctDoubleType_t(8, excit_type%fullstart_lowering, excit_type%fullstart_lowering, '1h'), & iikj = DistinctDoubleType_t(9, excit_type%fullstart_R_to_L, excit_type%non_overlap, '1i', [1, 4, 3, 2]), & iijk = DistinctDoubleType_t(10, excit_type%fullstart_L_to_R, excit_type%non_overlap, '1j', [1, 4, 3, 2]), & ijij = DistinctDoubleType_t(11, excit_type%fullstart_stop_alike, excit_type%fullstart_stop_alike, '2a'), & jiji = DistinctDoubleType_t(12, excit_type%fullstart_stop_alike, excit_type%fullstart_stop_alike, '2b'), & iijj = DistinctDoubleType_t(13, excit_type%fullstart_stop_mixed, excit_type%non_overlap, '2c', [1, 4, 3, 2]), & ikjl = DistinctDoubleType_t(14, excit_type%double_raising, excit_type%double_raising, '3a', [3, 4, 1, 2]), & kilj = DistinctDoubleType_t(15, excit_type%double_lowering, excit_type%double_lowering, '3b'), & ijkl = DistinctDoubleType_t(16, excit_type%double_R_to_L_to_R, excit_type%non_overlap, '3c', [1, 4, 3, 2]), & jilk = DistinctDoubleType_t(17, excit_type%double_L_to_R_to_L, excit_type%non_overlap, '3d', [1, 4, 3, 2]), & ijlk = DistinctDoubleType_t(18, excit_type%double_R_to_L, excit_type%non_overlap, '3e', [1, 4, 3, 2]), & jikl = DistinctDoubleType_t(19, excit_type%double_L_to_R, excit_type%non_overlap, '3f', [1, 4, 3, 2]) end type type(DistinctDoubleTypeVals_t), parameter :: distinct_doubles = DistinctDoubleTypeVals_t() type(DistinctDoubleType_t), parameter :: distinct_doubles_arr(0:19) & = [ & distinct_doubles%invalid, & distinct_doubles%jijk, & distinct_doubles%ijkj, & distinct_doubles%ikjk, & distinct_doubles%kikj, & distinct_doubles%jikk, & distinct_doubles%ijkk, & distinct_doubles%ijik, & distinct_doubles%jiki, & distinct_doubles%iikj, & distinct_doubles%iijk, & distinct_doubles%ijij, & distinct_doubles%jiji, & distinct_doubles%iijj, & distinct_doubles%ikjl, & distinct_doubles%kilj, & distinct_doubles%ijkl, & distinct_doubles%jilk, & distinct_doubles%ijlk, & distinct_doubles%jikl & ] interface DistinctDoubleType_t module procedure construct_DistinctDoubleType_t_ijkl, construct_DistinctDoubleType_t_ijab end interface type :: ijkl_Index_t !! An index in the format of the Shavit papers [@Hinze1981] \( e_{ijkl} \) !! !! This means j and l are the particles, i and k the holes. integer :: val(4) contains procedure :: extract_ijkl => extract_ijkl_Index_t procedure :: convert => convert_ijkl_Index_t procedure :: normalize => normalize_ijkl_Index_t procedure :: is_normalized => is_normalized_ijkl_Index_t procedure :: fuse => fuse_ijkl_Index_t end type type :: ijab_Index_t !! An index in the format of the rest of `NECI` !! [i, j, a, b] corresponds to \(e_{aibj} \) !! !! This means i, j are the particles, a and b the holes integer :: val(4) contains procedure :: extract_ijab => extract_ijab_Index_t procedure :: convert => convert_ijab_Index_t procedure :: normalize => normalize_ijab_Index_t procedure :: is_normalized => is_normalized_ijab_Index_t procedure :: fuse => fuse_ijab_Index_t end type type :: DistinctDouble_t type(ijkl_Index_t) :: idx type(DistinctDoubleType_t) :: typ contains private procedure, public :: get_repr_excit_info procedure, public :: get_repr_index procedure, public :: fuse => fuse_DistinctDouble_t end type interface DistinctDouble_t module procedure construct_DistinctDouble_t_ijkl, construct_DistinctDouble_t_ijab end interface type :: ExcitationInformation_t type(ExcitationType_t) :: typ = excit_type%invalid ! need the involved indices of the excitation: list of integers ! for now the convention is, that they are given in an ordered form ! and is not related to the involved generators E_{ij} (E_{kl}) ! directly, for single excitations ofc. only two entries of this ! vector needed. ! update: ! new convention store, original indiced and the ordered ones in ! the fullStart, etc. indices. integer :: i = -1 integer :: j = -1 integer :: k = -1 integer :: l = -1 integer :: fullStart = -1 integer :: fullEnd = -1 integer :: secondStart = -1 integer :: firstEnd = -1 integer :: weight = -1 ! can get rid of this in future! ! misuse secondstart firstend -> as weight as it is not used in the typ ! of excitations where weights is needed ! update: ! dont need overlaprange, and nonoverlaprange anywhere, just need to ! indicate if its a non-overlap, single overlap or proper double! integer :: overlap = -1 ! in rework get rid of this too! ! not needed, %typ could be used to get same results ! since eg. calcRemainingSwitches is only needed in cases, where there ! is atleast one overlap site or in single excitations! ! where there is no overlap at all! ! 0 ... no overlap or single ! 1 ... single overlap ! >1 ... proper double overlap ! generator flags: necessary information on involved generators. ! could store it as flags (0,1) but that would mask the meaning and ! also ruin some matrix accessing functionality: so for now store ! lowering: -1 ! raising : +1 ! weight: 0 ! maybe need firstGen, lastGen, (even secondGen) maybe? integer :: gen1 = -2 integer :: gen2 = -2 integer :: firstGen = -2 integer :: lastGen = -2 integer :: currentGen = -2 ! also store excitation level(number of occupation number differences) ! in this type, to have everything at one place integer :: excitLvl = -1 ! definetly get rid of that! never used ! at all! -> UPDATE! with changing of relative probabilities of ! those excitations, i definetly need this type of information! ! misuse it in such a way, that i store 5 different types of double ! excitations!: ! 0.. (ii,jj) RR/LL ! 1.. (ii,jj) RL ! 2.. (ii,jk) RR/LL ! 3.. (ii,jk) RL ! 4.. (ij,kl) x ! additional flags: ! for a 4 index double excitation, there is an additional flag ! necessary to distiguish between certain kind of equally possible ! excitations: ! logical :: fourFlag ! the order of generators is some excitations has an influence on the ! relative sign of the x_1 semi-stop matrix elements ! use a real here: 1.0 or -1.0 and just multiply x1 element real(dp) :: order = 0.0_dp real(dp) :: order1 = 0.0_dp ! maybe can get rid of order parameters.. since i could in general ! always choose such generators that the order parameter is +1 ! but that did not work beforhand... ? hm ! !TODO maybe more necessary info needed. ! add a flag to indicate if the excitation is even possible or other ! wise cancel the excitation logical :: valid = .false. ! for the exact calculation, to avoid calculating non-overlap ! contributions to matrix elements which are not possible, due to ! spin-coupling changes in the overlap range use a flag to ! indicate if a spin_change happened logical :: spin_change = .false. integer :: i_sg_start = -1 !! The supergroup from which the excitation starts. !! Only relevant for PropVec. end type ExcitationInformation_t ! logical to indicate that GUGA and core space, like doubles and singles ! are used as the semi-stochastic core space logical :: tGUGACore = .false. ! also use a memory-tag for the exact excitations array to keep track ! how big those arrays get.. and also of the bigger temp_excits integer(TagIntType) :: tag_excitations = 0, tag_tmp_excits = 0 ! also store the projected energy lists.. but only for 1 run not for ! all of them i guess. integer(TagIntType) :: tag_proje_list = 0 ! ======================== end TYPE defs ================================ ! ========================= INTERFACES ================================== ! create derived type to use array of procedure pointers to efficiently ! access the functions necessary for the matrix element calculation abstract interface pure function dummyFunction(bValue) result(ret) import :: dp implicit none real(dp), intent(in) :: bValue real(dp) :: ret end function dummyFunction pure subroutine dummySubroutine(bValue, x0, x1) import :: dp implicit none real(dp), intent(in) :: bValue real(dp), intent(out), optional :: x0, x1 end subroutine dummySubroutine end interface type :: ProcedurePtrArray_t procedure(dummyFunction), pointer, nopass :: ptr => null() ! why the nopass flag is needed see: ! https://software.intel.com/en-us/forums/topic/508530 end type ProcedurePtrArray_t type :: ProcPtrArrTwo_t procedure(dummySubroutine), pointer, nopass :: ptr => null() end type ProcPtrArrTwo_t ! all procedure pointers in an array have to use the same interface! ! no generic or differing ones for each element are allowed! -> this ! means that one has to use the most general one(3 input, 1 output) for ! all, even if everytime a 1 is given as output, as in many cases -> ! ask Simon if there is room for improvement here ! create an array of procedure pointers to all necessary functions for ! single excitation matrix element calculation type(ProcedurePtrArray_t) :: singleMatElesGUGA(15) ! define an array to give the correct indices to access matrix element ! function for given stepvectors, delta b, and generator type ! delta b value = {-1,1} and generator flags = {-1...lowering,1...raising} ! fill up with correct values, -1 corresponds to an impossible combination ! and will cause an error when tried to use as an index for matrix element ! calculation! are combined to a single index ranging from -1:2 ! by accessing through: db + (G+1)/2 integer, dimension(0:3,0:3,-1:2) :: indArrOne = reshape( (/ & !& 00,10,20, 30,01,11, 21,31,02,12, 22,32, 03,13,23,33: 2, 2, 2, -1, 2,10, 14, 5, 2, -1, 3, 7, -1, 7, 5, 3, & ! DeltaB = -1 & L 2, 2, 2, -1, 2, 3, 12, 6, 2, -1,11, 8, -1, 5, 7, 3, & ! DeltaB = -1 & R 2, 2, 2, -1, 2, 3, -1, 5, 2, 13,10, 7, -1, 7, 5, 3, & ! DeltaB = +1 & L 2, 2, 2, -1, 2, 9, -1, 6, 2, 15, 3, 8, -1, 5, 7, 3 & ! DeltaB = +1 & R /), (/ 4, 4, 4 /)) ! use two different arrays of procedure pointers for the x=0 and x=1 ! double excitation matrix element ! x1 elements: type(ProcedurePtrArray_t) :: doubleMatEleX1GUGA(45) ! and make access index-matrix: ! probably an error here: element 45-> links to non-existent function ! funB(b,3,2) and may be in the wrong position too.. integer, dimension(0:3,0:3,-7:7) :: indArrTwoX1 = reshape( (/ & !& !00,10, 20, 30, 01, 11, 21, 31, 02, 12, 22, 32, 03, 13, 23, 33: 2, -1, 2, -1, -1, 34, 32, 18, -1, -1, 2, -1, -1, -1, -1, 2, &! db = -2 & LL 2, -1, 2, -1, 2, 39, 43, 23, -1, -1, 39, -1, -1, -1, 15, 2, &! db = -2 & RL 2, -1, -1, -1, 2, 2, 33, -1, -1, -1, 45, -1, -1, -1, 3, 2, &! db = -2 & RR -1, -1, -1, -1, 16, -1, -1, -1, 2, -1, -1, -1, -1, 3, 7, -1, &! db= -1 & LL 1, 38, 7, 1, 12, 12, 40, 17, 38, 38, 7, 23, 1, 4, 20, 1, &! db = -1 & RL + full-start! -1, 2, 21, -1, -1, -1, -1, 12, -1, -1, -1, 22, -1, -1, -1, -1, &! db = -1 & RR->maybe store info here 2, 21, 16, 1, -1, 31, 33, 28, -1, 32, 35, 26, 1, -1, -1, 2, &! db = 0 & LL 2, 7, 12, -1, 28, 30, 44, 10, 26, 42, 34, 6, -1, 20, 17, 2, &! db = 0 & RL 2, -1, -1, 1, 11, 31, 37, -1, 5, 29, 35, -1, 1, 7, 12, 2, &! db = 0 & RR -1, -1, -1, -1, 2, -1, -1, -1, 21, -1, -1, -1, -1, 12, 8, -1, &! db = +1 & LL -1, 12, 40, -1, 40, -1, -1, 25, 7, -1, -1, 20, -1, 17, 9, -1, &! db = +1 & RL -1, 16, 2, -1, -1, -1, -1, 24, -1, -1, -1, 7, -1, -1, -1, -1, &! db = +1 & RR 2, 2, -1, -1, -1, 2, -1, -1, -1, 33, 30, 14, -1, -1, -1, 2, &! db = +2 & LL 2, 2, -1, -1, -1, 41, -1, -1, 2, 43, 41, 25, -1, 19, -1, 2, &! db = +2 & RL 2, -1, -1, -1, -1, 36, -1, -1, 2, 32, 2, -1, -1, 8, -1, 2 &! db = +2 & RR /), (/4, 4, 15/)) ! x0 elemets: type(ProcedurePtrArray_t) :: doubleMatEleX0GUGA(17) ! build index matrix ! could make some matrix elements indpendent of deltaB value, by filling ! every entry with the same matrix element, but for now choose to only ! make certain starts possible... ! also could point some impossible stepvector combinations to the ! zero function, and not making the index -1 -> decide on that! ! i could do smth like abort the matrix element calculation if its a -1 ! index in the matrix, like abort the calcuation and output 0.. ! REMEMBER: x0 matrix element only needed when Delta b = 0 all over excitation... ! maybe room for improvement there... -> finish that but think about that! integer, dimension(0:3,0:3,-7:7) :: indArrTwoX0 = reshape( (/ & !& !00, 10,20, 30, 01,11,21, 31, 02, 12,22, 32, 03, 13, 23,33: 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, & ! db=-2 & LL -> this line always leads to 0 matEle... 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, & ! db=-2 & RL -> also this 0.. 1, -1, -1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1, & ! db=-2 & RR -> this line is always 0... 1, -1, -1, -1, 14, 1, 2, -1, 1, -1, 2, -1, -1, 1, 7, 2, & ! db=-1 & LL 1, 1, 7, 10, 7, 7, 1, 15, 1, 1, 7, 1, 11, 1, 17, 5, & ! db=-1 & RL -> also store full-start elements here 1, 1, 16, -1, -1, 2, 2, 7, -1, -1, 1, 1, -1, -1, -1, 2, & ! db=-1 & RR -> do include single overlap matrix elements! 2, 16, 14, 4, -1, 3, 1, 7, -1, 1, 3, 7, 4, -1, -1, 2, & ! db= 0 & LL 2, 6, 6, -1, 6, 2, 1, 13, 6, 1, 2, 12, -1, 16, 14, 2, & ! db= 0 & RL 2, -1, -1, 4, 13, 3, 1, -1, 12, 1, 3, -1, 4, 7, 7, 2, & ! db= 0 & RR 1, -1, -1, -1, 1, 2,-1, -1, 16, 2, 1, -1, -1, 7, 1, 2, & ! db=+1 & LL -1, 7, 1, 9, 1,-1,-1, 1, 7, -1,-1, 17, 8, 15, 1,-1, & ! db=+1 & RL 1, 14, 1, -1, -1, 1,-1, 1, -1, 2, 2, 7, -1, -1, -1, 2, & ! db=+1 & RR 1, 1, -1, -1, -1, 1,-1, -1, -1, 1, 1, 1, -1, -1, -1, 1, & ! db=+2 & LL -> always 0 1, 1, -1, -1, -1, 1,-1, -1, 1, 1, 1, 1, -1, 1, -1, 1, & ! db=+2 & RL 1, -1, -1, -1, -1, 1,-1, -1, 1, 1, 1, -1, -1, 1, -1, 1 & ! db=+2 & RR -> always 0 /),(/ 4, 4, 15 /)) ! to index third matrix dimension use: (db + 2)*3 + (G1 + G2)/2 + 2 ! where G = +1 for R and -1 for L ! unfortunately for mixed generators, full-stop matrix element is needed ! seperately, because otherwise there would be ambiguities with RL ! intermediate matrix elements ! although procedure pointer arrays too much for this case.. as it is ! pretty special type(ProcPtrArrTwo_t) :: mixedGenFullStopMatEle(5) ! find good indexing converting function to efficiently access those ! functions. As only some combinations are needed, but the restriction on ! the delta b values has to be included here too. ! could write function or do matrix again, but almost only -1 in matrix... integer, dimension(0:3, 0:3, -1:1) :: indArrEnd = reshape([ & -1, -1, -1, -1, -1, -1, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, & ! db = -2 1, -1, -1, -1, -1, 2, -1, -1, -1, -1, 3, -1, -1, -1, -1, 4, & ! db = 0 -1, -1, -1, -1, -1, -1, -1, -1, -1, 5, -1, -1, -1, -1, -1, -1 & ! db = +2 ], [4, 4, 3]) ! and create a indexing array to acces matrix elements through: ! doubleMatEle(indArrTwo(d',d',db,G1,G2,x0,x1,b)) ! think about a good way to combine possible ! also need a similar procedure pointer array for the necessary calculation ! of the r_k terms, coming from the two-particle contributions to single ! excitations. they behave very similar to the usual single excitation ! product terms and hence, it is probably convenient to calculate them ! in a similar fashion. although it might be a overkill, since only 2 ! stepvector combinations correspond to a b-dependent function and the ! rest are just constants... think ybout that efficiency and ask simon it ! this is too costly.. type(ProcedurePtrArray_t) :: doubleContribution(7) ! write similar matrix indication table ! delete out the start and end values, as they are never needed! ! and should not be accessed! integer, dimension(0:3,0:3,-1:2) :: indContr = reshape( [ & !&! 00, 10, 20, 30, 01, 11, 21, 31, 02, 12, 22, 32, 03, 13, 23, 33 1, -1, -1, -1, -1, 1, 7, -1, -1, -1, 3, -1, -1, -1, -1, 3, &! db = -1 & L 1, -1, -1, -1, -1, 3, 5, -1, -1, -1, 1, -1, -1, -1, -1, 3, &! db = -1 & R 1, -1, -1, -1, -1, 3, -1, -1, -1, 6, 1, -1, -1, -1, -1, 3, &! db = +1 & L 1, -1, -1, -1, -1, 1, -1, -1, -1, 4, 3, -1, -1, -1, -1, 3 &! db = +1 & R ], [4,4,4]) ! also need a list of determinants and matrix elements connceted to the ! reference determinant ! adapt that to multiple neci runs.. hope that works as intended.. ! probably have to use it as a type, to store lists or different lists ! in it, and also be able to (de)allocate them individually type ProjE_t integer(n_int), allocatable :: projE_ilut_list(:, :) HElement_t(dp), allocatable :: projE_hel_list(:) ! also store the excitation level in the projected list, since otherwise ! it is really hard to determine it in the GUGA formalism integer, allocatable :: exlevel(:) ! also store the number of entries to correctly binary search integer :: num_entries end type ProjE_t type(ProjE_t), allocatable :: projE_replica(:) type RdmContribList_t integer(n_int), allocatable :: ilut_list(:, :) real(dp), allocatable :: rdm_contrib(:) integer(int_rdm), allocatable :: rdm_ind(:) integer, allocatable :: repeat_count(:) end type RdmContribList_t ! also make a global integer list of orbital indices, so i do not have to ! remake them in every random orbital picker! integer, allocatable :: orbitalIndex(:) contains subroutine init_guga_data_procPtrs() ! this subroutine initializes the procedure pointer arrays needed for ! the matrix element calculation ! -------- single excitation procedure pointers:--------------------- ! and point them to the corresponding functions -> for now write down all ! differing functions explicetly, as i can't think of somw fancy way to ! include that better ! store only procedure pointer to differing functions and convert to ! indices in such a way that correct ones get picked! ! also store weight-generator matrix element, as it may be of some use singleMatElesGUGA(1)%ptr => funZero singleMatElesGUGA(2)%ptr => funPlus1 singleMatElesGUGA(3)%ptr => funMinus1 singleMatElesGUGA(4)%ptr => funTwo singleMatElesGUGA(5)%ptr => funA_0_1 singleMatElesGUGA(6)%ptr => funA_1_0 singleMatElesGUGA(7)%ptr => funA_2_1 singleMatElesGUGA(8)%ptr => funA_1_2 singleMatElesGUGA(9)%ptr => funC_0 singleMatElesGUGA(10)%ptr => funC_1 singleMatElesGUGA(11)%ptr => funC_2 singleMatElesGUGA(12)%ptr => funOverB_0 singleMatElesGUGA(13)%ptr => funOverB_1 singleMatElesGUGA(14)%ptr => minFunOverB_1 singleMatElesGUGA(15)%ptr => minFunOverB_2 ! ----------- double excitation procedure pointer -------------------- ! ------ double excitation x0 matrix element part -------------------- ! and fill it up with all unique values needed for x=1 matrix elements doubleMatEleX1GUGA(1)%ptr => funZero doubleMatEleX1GUGA(2)%ptr => funPlus1 ! doubleMatEleX1GUGA(3)%ptr => funMinus1 doubleMatEleX1GUGA(3)%ptr => funA_3_2 doubleMatEleX1GUGA(4)%ptr => minFunA_3_2 doubleMatEleX1GUGA(5)%ptr => funA_3_2overR2 doubleMatEleX1GUGA(6)%ptr => minFunA_3_2overR2 ! doubleMatEleX1GUGA(8)%ptr => funA_0_2overR2 doubleMatEleX1GUGA(7)%ptr => minFunA_0_2_overR2 doubleMatEleX1GUGA(8)%ptr => funA_m1_0 doubleMatEleX1GUGA(9)%ptr => minFunA_m1_0 doubleMatEleX1GUGA(10)%ptr => funA_m1_0_overR2 doubleMatEleX1GUGA(11)%ptr => minFunA_m1_0_overR2 doubleMatEleX1GUGA(12)%ptr => funA_2_0_overR2 doubleMatEleX1GUGA(13)%ptr => minFunA_2_0_overR2 doubleMatEleX1GUGA(14)%ptr => funA_2_1 doubleMatEleX1GUGA(15)%ptr => minFunA_2_1 doubleMatEleX1GUGA(16)%ptr => funA_2_1_overR2 doubleMatEleX1GUGA(17)%ptr => minFunA_2_1_overR2 doubleMatEleX1GUGA(18)%ptr => funA_0_1 doubleMatEleX1GUGA(19)%ptr => minFunA_0_1 doubleMatEleX1GUGA(20)%ptr => funA_0_1_overR2 doubleMatEleX1GUGA(21)%ptr => minFunA_0_1_overR2 doubleMatEleX1GUGA(22)%ptr => funA_1_2 doubleMatEleX1GUGA(23)%ptr => minFunA_1_2 doubleMatEleX1GUGA(24)%ptr => funA_1_0 doubleMatEleX1GUGA(25)%ptr => minFunA_1_0 doubleMatEleX1GUGA(26)%ptr => funA_3_1_overR2 doubleMatEleX1GUGA(27)%ptr => minFunA_3_1_overR2 ! doubleMatEleX1GUGA(30)%ptr => funA_m1_1_overR2 doubleMatEleX1GUGA(28)%ptr => minFunA_m1_1_overR2 doubleMatEleX1GUGA(29)%ptr => funB_2_3 doubleMatEleX1GUGA(30)%ptr => funD_0 doubleMatEleX1GUGA(31)%ptr => minFunD_0 doubleMatEleX1GUGA(32)%ptr => funB_1_2 doubleMatEleX1GUGA(33)%ptr => funB_0_1 doubleMatEleX1GUGA(34)%ptr => funD_1 doubleMatEleX1GUGA(35)%ptr => minFunD_1 doubleMatEleX1GUGA(36)%ptr => funD_m1 doubleMatEleX1GUGA(37)%ptr => funB_m1_0 doubleMatEleX1GUGA(38)%ptr => funC_2 doubleMatEleX1GUGA(39)%ptr => minFunC_2 doubleMatEleX1GUGA(40)%ptr => funC_0 doubleMatEleX1GUGA(41)%ptr => minFunC_0 doubleMatEleX1GUGA(42)%ptr => minFunOverB_2_R2 doubleMatEleX1GUGA(43)%ptr => minFunB_0_2 doubleMatEleX1GUGA(44)%ptr => minFunOverB_0_R2 doubleMatEleX1GUGA(45)%ptr => funD_2 ! -------------- double excitation x0 matrix element part ------------ doubleMatEleX0GUGA(1)%ptr => funZero doubleMatEleX0GUGA(2)%ptr => funPlus1 doubleMatEleX0GUGA(3)%ptr => funMinus1 doubleMatEleX0GUGA(4)%ptr => funSqrt2 doubleMatEleX0GUGA(5)%ptr => minFunSqrt2 doubleMatEleX0GUGA(6)%ptr => funOverRoot2 doubleMatEleX0GUGA(7)%ptr => minFunOverR2 doubleMatEleX0GUGA(8)%ptr => funA_0_1 doubleMatEleX0GUGA(9)%ptr => funA_1_0 doubleMatEleX0GUGA(10)%ptr => funA_1_2 doubleMatEleX0GUGA(11)%ptr => funA_2_1 doubleMatEleX0GUGA(12)%ptr => funA_1_2overR2 doubleMatEleX0GUGA(13)%ptr => funA_1_0_overR2 doubleMatEleX0GUGA(14)%ptr => funA_0_1_overR2 doubleMatEleX0GUGA(15)%ptr => minFunA_0_1_overR2 doubleMatEleX0GUGA(16)%ptr => funA_2_1_overR2 doubleMatEleX0GUGA(17)%ptr => minFunA_2_1_overR2 ! ------ mixed generator full-stop matrix elements-------------------- mixedGenFullStopMatEle(1)%ptr => fullStop_00 mixedGenFullStopMatEle(2)%ptr => fullStop_11 mixedGenFullStopMatEle(3)%ptr => fullStop_22 mixedGenFullStopMatEle(4)%ptr => fullStop_33 mixedGenFullStopMatEle(5)%ptr => fullStop_12 ! --------- double contributions to single excitaitons---------------- doubleContribution(1)%ptr => funZero doubleContribution(2)%ptr => funPlus1 doubleContribution(3)%ptr => funMinus1 doubleContribution(4)%ptr => minFunBplus2 doubleContribution(5)%ptr => funBplus0 doubleContribution(6)%ptr => funBplus1 doubleContribution(7)%ptr => minFunBplus1 end subroutine init_guga_data_procPtrs ! wrapper functions to access matrix element terms ! for consistency reasons also write an get.. function for the specific ! mixed gen full-stops. elemental subroutine getMixedFullStop(step1, step2, deltaB, bValue, x0_element, & x1_element) ! function to access the special mixed generator full-stop elements ! which due to storage reasons are stored in a seperate func. pointer ! array integer, intent(in) :: step1, step2, deltaB real(dp), intent(in) :: bValue real(dp), intent(out), optional :: x0_element, x1_element integer :: ind ! get index: ind = indArrEnd(step1, step2, deltaB / 2) ! with the optional output arguments can also just calc. x0 or x1 call mixedGenFullStopMatEle(ind)%ptr(bValue, x0_element, x1_element) end subroutine getMixedFullStop elemental function getDoubleContribution(step1, step2, deltaB, genFlag, bValue) & result(doubleContr) ! Access necessary two-particle contribution to single excitation ! matrix elements. ! ! input: ! step1/2 ... stepvector values of CSFs, step1 is from <d'| ! deltaB ... current delta b value of excitation ! genFlag ... generator flag of excitation ! bValue ... current b value of CSF ! ! output: ! doubleContr ... product term from shavitt graph rules used to calculate ! the matrix element between two given CSFs integer, intent(in) :: step1, step2, deltaB, genFlag real(dp), intent(in) :: bValue real(dp) :: doubleContr integer :: ind ! get index ind = indContr(step1, step2, deltaB + (genFlag + 1) / 2) doubleContr = doubleContribution(ind)%ptr(bValue) end function getDoubleContribution elemental function getSingleMatrixElement(step1, step2, deltaB, genFlag, bValue) & result(hElement) ! Access the necessary single excitation product terms for the H matrix ! element calculation. ! ! input: ! step1/2 ... stepvector values of CSFs, step1 is from <d'| ! deltaB ... current delta b value of excitation ! genFlag ... generator flag of excitation ! bValue ... current b value of CSF ! ! output: ! hElement ... product term from shavitt graph rules used to calculate ! the matrix element between two given CSFs integer, intent(in) :: step1, step2, deltaB, genFlag real(dp), intent(in) :: bValue real(dp) :: hElement integer :: ind ! only need for this function is to correctly access the procedure ! pointers to get correct function -> get index from index matrix: ind = indArrOne(step1, step2, deltaB + (genFlag + 1) / 2) ! call correct function: hElement = singleMatElesGUGA(ind)%ptr(bValue) end function getSingleMatrixElement elemental subroutine getDoubleMatrixElement(step1, step2, deltaB, genFlag1, genFlag2, & bValue, order, x0_element, x1_element) ! access the necessary double excitation product terms for the H matrix ! element calculation ! ! input: ! step1/2 ... stepvector values if given CSFs, step1 form <d'| ! deltaB ... current delta b value of excitation ! genFlag1/2... generator types involved in excitation ! bValue ... current b value of CSF ! order ... order parameter determining the sign of x1 element ! ! output: ! x0_element... x0 product term (optional) ! x1_element... x1 product term (optional) ! -> call this function by x0_element=var, x1_element=var2, if only ! specific matrix element is wanted. x0_element is default if only ! one output variable is present integer, intent(in) :: step1, step2, deltaB, genFlag1, genFlag2 real(dp), intent(in) :: bValue, order real(dp), intent(out), optional :: x0_element, x1_element integer :: x0_ind, x1_ind if (present(x0_element)) then ! first get correct indices to access procedure pointer array x0_ind = indArrTwoX0(step1, step2, 3 * deltaB + (genFlag1 + genFlag2) / 2) ! then call corresponding function x0_element = doubleMatEleX0GUGA(x0_ind)%ptr(bValue) end if ! same for x1 element if (present(x1_element)) then x1_ind = indArrTwoX1(step1, step2, 3 * deltaB + (genFlag1 + genFlag2) / 2) x1_element = order * doubleMatEleX1GUGA(x1_ind)%ptr(bValue) end if end subroutine getDoubleMatrixElement ! maybe write specific functions to only access x0/x1 matrix elements to ! avoid double conditioning... ! ===== special functions for mixed generator full-stop elements ======== ! maybe create subroutines to output x0 and x1 matrix elements for all ! needed terms, and let the procedure pointers point to them pure subroutine fullStop_00(b, x0, x1) real(dp), intent(in) :: b real(dp), intent(out), optional :: x0, x1 unused_var(b) if (present(x0)) x0 = 0.0_dp if (present(x1)) x1 = 0.0_dp end subroutine fullStop_00 pure subroutine fullStop_11(b, x0, x1) real(dp), intent(in) :: b real(dp), intent(out), optional :: x0, x1 if (present(x0)) x0 = OverR2 if (present(x1)) x1 = minFunA_m1_1_overR2(b) end subroutine fullStop_11 pure subroutine fullStop_12(b, x0, x1) ! is the same for switched stepvector values real(dp), intent(in) :: b real(dp), intent(out), optional :: x0, x1 unused_var(b) if (present(x0)) x0 = 0.0_dp if (present(x1)) x1 = 1.0_dp end subroutine fullStop_12 pure subroutine fullStop_22(b, x0, x1) real(dp), intent(in) :: b real(dp), intent(out), optional :: x0, x1 if (present(x0)) x0 = OverR2 if (present(x1)) x1 = funA_3_1_overR2(b) end subroutine fullStop_22 pure subroutine fullStop_33(b, x0, x1) real(dp), intent(in) :: b real(dp), intent(out), optional :: x0, x1 unused_var(b) if (present(x0)) x0 = Root2 if (present(x1)) x1 = 0.0_dp end subroutine fullStop_33 ! ===== special functions for the double contribution to single excitation ! matrix elements pure function minFunBplus2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -(b + 2.0_dp) end function minFunBplus2 pure function funBplus0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = b end function funBplus0 pure function minFunBplus1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -(b + 1.0_dp) end function minFunBplus1 pure function funBplus1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = b + 1.0_dp end function funBplus1 ! =========== additional double excitation matrix elements =============== pure function funSqrt2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = Root2 end function funSqrt2 pure function minFunSqrt2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = -Root2 end function minFunSqrt2 pure function funOverRoot2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = OverR2 end function funOverRoot2 pure function minFunOverR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = -OverR2 end function minFunOverR2 pure function funA_1_2overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA_1_2(b) end function funA_1_2overR2 pure function funA_3_2overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA_3_2(b) end function funA_3_2overR2 pure function minFunA_3_2overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_3_2overR2(b) end function minFunA_3_2overR2 pure function funA_0_2overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA(b, 0.0_dp, 2.0_dp) end function funA_0_2overR2 pure function minFunA_0_2_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_0_2overR2(b) end function minFunA_0_2_overR2 pure function funA_3_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funA(b, 3.0_dp, 2.0_dp) end function funA_3_2 pure function minFunA_3_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_3_2(b) end function minFunA_3_2 pure function funA_1_0_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA_1_0(b) end function funA_1_0_overR2 pure function funA_m1_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funA(b, -1.0_dp, 0.0_dp) end function funA_m1_0 pure function minFunA_m1_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_m1_0(b) end function minFunA_m1_0 pure function funA_m1_0_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA_m1_0(b) end function funA_m1_0_overR2 pure function minFunA_m1_0_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_m1_0_overR2(b) end function minFunA_m1_0_overR2 pure function funA_2_0_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA(b, 2.0_dp, 0.0_dp) end function funA_2_0_overR2 pure function minFunA_2_0_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_2_0_overR2(b) end function minFunA_2_0_overR2 pure function funA_0_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA_0_1(b) end function funA_0_1_overR2 pure function minFunA_0_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_0_1_overR2(b) end function minFunA_0_1_overR2 pure function funA_2_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA_2_1(b) end function funA_2_1_overR2 pure function minFunA_2_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_2_1_overR2(b) end function minFunA_2_1_overR2 pure function minFunA_1_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_1_2(b) end function minFunA_1_2 pure function minFunA_1_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_1_0(b) end function minFunA_1_0 pure function minFunA_0_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_0_1(b) end function minFunA_0_1 pure function funA_3_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA(b, 3.0_dp, 1.0_dp) end function funA_3_1_overR2 pure function minFunA_3_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_3_1_overR2(b) end function minFunA_3_1_overR2 pure function funA_m1_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = OverR2 * funA(b, -1.0_dp, 1.0_dp) end function funA_m1_1_overR2 pure function minFunA_m1_1_overR2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_m1_1_overR2(b) end function minFunA_m1_1_overR2 pure function funB_2_3(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funB(b, 2.0_dp, 3.0_dp) end function funB_2_3 pure function funD_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funD(b, 2.0_dp) end function funD_2 pure function funD_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funD(b, 0.0_dp) end function funD_0 pure function minFunD_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funD_0(b) end function minFunD_0 pure function funB_1_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funB(b, 1.0_dp, 2.0_dp) end function funB_1_2 pure function funB_0_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funB(b, 0.0_dp, 1.0_dp) end function funB_0_1 pure function funD_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funD(b, 1.0_dp) end function funD_1 pure function minFunD_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funD_1(b) end function minFunD_1 pure function funD_m1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funD(b, -1.0_dp) end function funD_m1 pure function funB_m1_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funB(b, -1.0_dp, 0.0_dp) end function funB_m1_0 ! mixed generator additional functions: pure function minFunA_2_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funA_2_1(b) end function minFunA_2_1 pure function minFunC_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funC_2(b) end function minFunC_2 pure function minFunC_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funC_0(b) end function minFunC_0 pure function minFunOverB_2_R2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = Root2 * minFunOverB_2(b) end function minFunOverB_2_R2 pure function minFunOverB_0_R2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -Root2 * funOverB_0(b) end function minFunOverB_0_R2 pure function minFunB_0_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funB(b, 0.0_dp, 2.0_dp) end function minFunB_0_2 !========= function for the single particle matrix calculation =========== ! ASSERT() probably not usable in "elemental" function, due to side-effects pure function funZero(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = 0.0_dp end function funZero pure function funPlus1(b) result(ret) ! think of a better way to include that -> wasted time real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = 1.0_dp end function funPlus1 pure function funMinus1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = -1.0_dp end function funMinus1 pure function funTwo(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret unused_var(b) ret = 2.0_dp end function funTwo pure function funA_0_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funA(b, 0.0_dp, 1.0_dp) end function funA_0_1 pure function funA_2_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funA(b, 2.0_dp, 1.0_dp) end function funA_2_1 pure function funA_1_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funA(b, 1.0_dp, 0.0_dp) end function funA_1_0 pure function funA_1_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funA(b, 1.0_dp, 2.0_dp) end function funA_1_2 pure function funC_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funC(b, 0.0_dp) end function funC_0 pure function funC_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funC(b, 1.0_dp) end function funC_1 pure function funC_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funC(b, 2.0_dp) end function funC_2 pure function minFunOverB_2(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funOverB(b, 2.0_dp) end function minFunOverB_2 pure function funOverB_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = funOverB(b, 1.0_dp) end function funOverB_1 pure function minFunOverB_1(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret ret = -funOverB(b, 1.0_dp) end function minFunOverB_1 pure function funOverB_0(b) result(ret) real(dp), intent(in) :: b real(dp) :: ret !ASSERT( b > 0.0_dp) ret = funOverB(b, 0.0_dp) end function funOverB_0 ! ==== end of specific single particle matrix element terms ============= ! =========== generic necessary functions =============================== ! use of ASSERT() probably not possible here, as it causes side-effects! pure function funA(b, x, y) result(ret) real(dp), intent(in) :: b, x, y real(dp) :: ret debug_function_name("funA") ASSERT( (b + y) > 0.0_dp) ret = sqrt((b + x) / (b + y)) end function funA pure function funB(b, x, y) result(ret) real(dp), intent(in) :: b, x, y real(dp) :: ret debug_function_name("funB") ASSERT( (b + x) > 0.0_dp) ASSERT( (b + y) > 0.0_dp) ret = sqrt(2.0_dp / ((b + x) * (b + y))) end function funB pure function funC(b, x) result(ret) real(dp), intent(in) :: b, x real(dp) :: ret debug_function_name("funC") ASSERT( (b + x) > 0.0_dp) ret = funA(b, x - 1.0_dp, x) * funA(b, x + 1.0_dp, x) end function funC pure function funD(b, x) result(ret) real(dp), intent(in) :: b, x real(dp) :: ret debug_function_name("funD") ASSERT( (b + x) > 0.0_dp) ret = funA(b, x + 2.0_dp, x) * funA(b, x - 1.0_dp, x + 1.0_dp) end function funD pure function funOverB(b, x) result(ret) real(dp), intent(in) :: b, x real(dp) :: ret debug_function_name("funOverB") ASSERT( (b + x) > 0.0_dp) ret = 1.0_dp / (b + x) end function funOverB ! =========== end of generic necessary functions ======================== elemental function construct_ExcitationType_t(n) result(res) integer, intent(in) :: n type(ExcitationType_t) :: res routine_name("construct_ExcitationType_t") res = excit_type_arr(n) ASSERT(res%val == n) end function elemental function get_repr_index(this) result(res) !! Return the index for the representation. !! !! If we have for example the `DistinctDoubleType_t` `jilk`, !! then the index of the representation shall refer !! to the `double_L_to_R_to_L` excitation and be `jkli` !! instead of the `non_overlap` part which would be `jilk`. class(DistinctDouble_t), intent(in) :: this type(ijkl_Index_t) :: res res = ijkl_Index_t(this%idx%val(this%typ%idx_for_repr)) end function elemental function get_repr_excit_info(this, i_sg_start) result(res) !! Return the index for the representation. !! !! If we have for example the `DistinctDoubleType_t` `jilk`, !! then the index of the representation shall refer !! to the `double_L_to_R_to_L` excitation and be `jkli` !! instead of the `non_overlap` part which would be `jilk`. class(DistinctDouble_t), intent(in) :: this integer, intent(in) :: i_sg_start type(ExcitationInformation_t) :: res type(ijkl_Index_t) :: idx integer :: i, j, k, l idx = this%get_repr_index() call idx%extract_ijkl(i, j, k, l) res = ExcitationInformation_t(this%typ%repr, i, j, k, l, i_sg_start=i_sg_start) end function elemental function construct_DistinctDoubleType_t_ijkl(idx) result(res) type(ijkl_Index_t), intent(in) :: idx type(DistinctDoubleType_t) :: res routine_name("construct_DistinctDoubleType_t_ijkl") integer :: i, j, k, l ASSERT(idx%is_normalized()) call idx%extract_ijkl(i, j, k, l) associate(dd => distinct_doubles) res = dd%invalid if (i == j) then if (i /= k .and. i /= l) then if (k == l) then res = dd%iijj else if (k < l) then res = dd%iijk else if (k > l) then res = dd%iikj end if end if else if (i < j) then if (i == k) then if (j == l) then res = dd%ijij else ASSERT(j < l) res = dd%ijik end if else if (i < k) then if (j < k) then if (j == l) then res = dd%ijkj else if (k == l) then res = dd%ijkk else if (k < l) then res = dd%ijkl else if (k > l) then res = dd%ijlk end if else if (j /= k) then if (j == l) then res = dd%ikjk else res = dd%ikjl end if end if end if end if else ASSERT (i > j) if (i == k) then if (j == l) then res = dd%jiji else if (k < l) then res = dd%jijk else if (k > l) then res = dd%kikj end if else if (i < k) then if (j == l) then res = dd%jiki else if (k == l) then res = dd%jikk else if (k < l) then res = dd%jikl else if (k > l) then if (i < l) then res = dd%jilk else if (i > l) then res = dd%kilj end if end if end if end if end associate end function elemental function construct_DistinctDoubleType_t_ijab(idx) result(res) type(ijab_Index_t), intent(in) :: idx type(DistinctDoubleType_t) :: res res = DistinctDoubleType_t(idx%convert()) end function elemental function construct_DistinctDouble_t_ijkl(idx) result(res) type(ijkl_Index_t), intent(in) :: idx type(DistinctDouble_t) :: res debug_function_name("construct_DistinctDouble_t_ijkl") ASSERT(idx%is_normalized()) res = DistinctDouble_t(idx, DistinctDoubleType_t(idx)) end function elemental function construct_DistinctDouble_t_ijab(idx) result(res) type(ijab_Index_t), intent(in) :: idx type(DistinctDouble_t) :: res debug_function_name("construct_DistinctDouble_t_ijkl") ASSERT(idx%is_normalized()) res = DistinctDouble_t(idx%convert()) end function elemental function normalize_ijab_Index_t(this) result(res) !! Normalize such that `i <= j` and `a <= b` class(ijab_Index_t), intent(in) :: this type(ijab_Index_t) :: res if (this%val(1) <= this%val(2)) then res%val(1:2) = this%val(1:2) else res%val(1:2) = this%val(2:1:-1) end if if (this%val(3) <= this%val(4)) then res%val(3:4) = this%val(3:4) else res%val(3:4) = this%val(4:3:-1) end if end function elemental function normalize_ijkl_Index_t(this) result(res) !! Normalize such that `i <= k` and `j <= l` class(ijkl_Index_t), intent(in) :: this type(ijkl_Index_t) :: res if (this%val(1) <= this%val(3)) then res%val(1) = this%val(1) res%val(3) = this%val(3) else res%val(1) = this%val(3) res%val(3) = this%val(1) end if if (this%val(2) <= this%val(4)) then res%val(2) = this%val(2) res%val(4) = this%val(4) else res%val(2) = this%val(4) res%val(4) = this%val(2) end if end function elemental function convert_ijkl_Index_t(this) result(res) !! Convert to the ijab_Index_t class(ijkl_Index_t), intent(in) :: this type(ijab_Index_t) :: res res = ijab_Index_t([this%val(2), this%val(4), this%val(1), this%val(3)]) end function elemental function convert_ijab_Index_t(this) result(res) !! Convert to the ijkl_Index_t class(ijab_Index_t), intent(in) :: this type(ijkl_Index_t) :: res res = ijkl_Index_t([this%val(3), this%val(1), this%val(4), this%val(2)]) end function elemental subroutine extract_ijab_Index_t(this, i, j, a, b) !! Extract the values as integers class(ijab_Index_t), intent(in) :: this integer, intent(out) :: i, j, a, b i = this%val(1) j = this%val(2) a = this%val(3) b = this%val(4) end subroutine elemental subroutine extract_ijkl_Index_t(this, i, j, k, l) !! Extract the values as integers class(ijkl_Index_t), intent(in) :: this integer, intent(out) :: i, j, k, l i = this%val(1) j = this%val(2) k = this%val(3) l = this%val(4) end subroutine elemental function is_normalized_ijab_Index_t(this) result(res) !! Check if index is normalized class(ijab_Index_t), intent(in) :: this logical :: res res = this%val(1) <= this%val(2) .and. this%val(3) <= this%val(4) end function elemental function is_normalized_ijkl_Index_t(this) result(res) !! Check if index is normalized class(ijkl_Index_t), intent(in) :: this logical :: res res = this%val(1) <= this%val(3) .and. this%val(2) <= this%val(4) end function elemental function fuse_ijkl_Index_t(this) result(res) !! Fuse a normalized index `i <= k` and `j <= l`. class(ijkl_Index_t), intent(in) :: this integer(int64) :: res associate(i => int(this%val(1), int64), j => int(this%val(2), int64), & k => int(this%val(3), int64), l => int(this%val(4), int64)) res = fuse_idx(fuse_symm_idx(i, k), fuse_symm_idx(j, l), ij_pairs, fortran=.false.) end associate end function elemental function fuse_ijab_Index_t(this) result(res) !! Fuse a normalized index `i <= j` and `a <= b`. class(ijab_Index_t), intent(in) :: this integer(int64) :: res associate(i => int(this%val(1), int64), j => int(this%val(2), int64), & a => int(this%val(3), int64), b => int(this%val(4), int64)) res = fuse_idx(fuse_symm_idx(a, b), fuse_symm_idx(i, j), ij_pairs, fortran=.false.) end associate end function elemental function fused_to_ijkl(ijkl) result(idx) !! Return a normalized index from a fused integer index integer(int64), intent(in) :: ijkl type(ijkl_Index_t) :: idx integer(int64):: i, j, k, l, ik, jl debug_function_name("fused_to_ijkl") call inv_fuse_idx(ijkl, ij_pairs, ik, jl, fortran=.false.) call inv_fuse_symm_idx(ik, i, k) call inv_fuse_symm_idx(jl, j, l) idx = ijkl_Index_t(int([i, j, k, l])) ASSERT(idx%is_normalized()) end function elemental function fused_to_ijab(ijab) result(idx) !! Return a normalized index from a fused integer index integer(int64), intent(in) :: ijab type(ijab_Index_t) :: idx integer(int64):: i, j, a, b, ij, ab debug_function_name("fused_to_ijab") call inv_fuse_idx(ijab, ij_pairs, ab, ij, fortran=.false.) call inv_fuse_symm_idx(ab, a, b) call inv_fuse_symm_idx(ij, i, j) idx = ijab_Index_t(int([i, j, a, b])) ASSERT(idx%is_normalized()) end function elemental function fuse_DistinctDouble_t(this) result(res) !! Fuse a Distinct Double excitation !! !! This means that a composite index for !! [typ, fused_ijkl_index] !! is created. !! !! @note !! This function does not preserve lexicographical ordering. !! I.e. if `(i_1, j_1, k_1, l_1) <= (i_2, j_2, k_2, l_2)` !! then the fused index is not necessarily `<=`. !! @endnote class(DistinctDouble_t), intent(in) :: this integer(int64) :: res res = fuse_idx(this%idx%fuse(), int(this%typ%val, int64) + 1, & size(distinct_doubles_arr, kind=int64), fortran=.false.) end function elemental function fused_to_DistinctDouble(n) result(distinct_double) !! Return a Distinct Double excitation from a fused integer index !! !! @note !! Note the set of valid `n` is not contiguos. !! @endnote integer(int64), intent(in) :: n type(DistinctDouble_t) :: distinct_double integer(int64) :: idx_typ, ijkl debug_function_name("fused_to_DistinctDouble") call inv_fuse_idx(n, size(distinct_doubles_arr, kind=int64), ijkl, idx_typ, fortran=.false.) distinct_double = DistinctDouble_t(& fused_to_ijkl(ijkl), & distinct_doubles_arr(idx_typ - 1) & ) #ifdef DEBUG_ associate(dd => DistinctDouble_t(fused_to_ijkl(ijkl))) ASSERT(dd%fuse() == n) end associate #endif end function end module