Allow specifying spec. orbitals for Breit

src/breit/CMakeLists.txt

@@ -24,4 +24,4 @@ add_library(rciqed-breit STATIC
     zkf.f90 zkf_I.f90
 )
 setup_fortran_modules(rciqed-breit)
-target_link_libraries_Fortran(rciqed-breit PRIVATE mod 9290 rang90)
+target_link_libraries_Fortran(rciqed-breit PRIVATE mod 9290 rang90 rciqed-grasp)

src/breit/bessel.f90

@@ -51,7 +51,7 @@ SUBROUTINE BESSEL(IA, IB, IK, IW, K)
       USE grid_C
       USE orb_C
       USE stor_C
-      USE wfaC_c
+      use grasp_rciqed_breit, only: WFACT, breit_specorbs
       IMPLICIT NONE
 !-----------------------------------------------
 !   D u m m y   A r g u m e n t s
@@ -85,7 +85,11 @@ SUBROUTINE BESSEL(IA, IB, IK, IW, K)
 !
 !   Function not in position; is it available in BESS arrays?
 !
-      W = WFACT*DABS(E(IA)-E(IB))/C
+      W = DABS(E(IA)-E(IB))/C
+      ! If either of the orbitals is designate as "non-spectroscopic" for the breit purposes,
+      ! we'll dampen the energy difference with WFACT (a very small factor) to effective go
+      ! to the frequency-independent limit of Breit
+      if(.not.breit_specorbs(IA) .or. .not.breit_specorbs(IB)) W = WFACT*W
       WIJ(IW) = W
 !
       DO IWKP = 1, 2

src/breit/breit.f90

@@ -1,40 +1,85 @@
-!> Routines related to calculating the matrix elements of the Breit operator.
+!> Routines related to calculating the matrix elements of the transverse photon operator
+!! (in the Coulomb gauge):
+!!
+!! \f[
+!!   g^T(R; \omega) =
+!!   - \vec{\alpha}_1 \cdot \vec{\alpha}_2 ~ \frac{\textrm{e}^{i\omega R}}{R}
+!!   - (\vec{\alpha}_1 \cdot \nabla_R) (\vec{\alpha}_2 \cdot \nabla_R) \frac{\textrm{e}^{i\omega R} - 1}{\omega^2 R}
+!! \f]
+!!
+!! Taking the frequency-independent limit \f$\lim_{\omega \to 0}\f$ leads to the Breit
+!! operator:
+!!
+!! \f[
+!!   g^B(R) = \lim_{\omega \to 0} g^T(R; \omega) =
+!!   - \frac{1}{2R} \left[
+!!     \vec{\alpha}_1 \cdot \vec{\alpha}_2
+!!   + (\vec{\alpha}_1 \cdot \hat{\vec{R}}) (\vec{\alpha}_2 \cdot \hat{\vec{R}})
+!!   \right]
+!! \f]
+!!
+!! In GRASP, this is not implemented directory, but with a \f$10^{-6}\f$ energy scaling factor
+!! `WFACT` (i.e. the value for \f$\omega\f$, derived from the difference of orbital energies,
+!! is multiplied by this factor).
+!!
+!! For completeness, we'll also mention another related operator: the Gaunt operator:
+!!
+!! \f[
+!!   g^G(R) = - \frac{\vec{\alpha}_1 \cdot \vec{\alpha}_2}{R}
+!! \f]
+!!
+!! While not implemented separately in GRASP, it is a different limit of the transverse
+!! photon operator.
 module grasp_rciqed_breit
     use, intrinsic :: iso_fortran_env, only: real64, dp => real64
+    use parameter_def, only: NNNW
     implicit none
 
+    !> Factor used to multiply the energy differences in the Breit operator to emulate
+    !! frequency-independent Breit.
+    real(real64), parameter :: WFACT = 1e-6_dp
+
+    !> Marks whether an orbital is considered to be a spectroscopic orbital or not for the
+    !! purposes of the Breit operator.
+    !!
+    !! The Breit two-particle matrix elements that contain at least one non-spectroscopic
+    !! orbital have their energy difference \f$\omega\f$ damped by `WFACT`.
+    logical :: breit_specorbs(NNNW)
+
 contains
 
     !> Initialize Breit-related global state.
     !!
     !! @param j2max This is the value determined by genintrk(), e.g. in init_rkintc()
-    subroutine init_breit(j2max)
+    !! @param specorbs List of `logical`s marking which orbitals are to be considered
+    !!        spectroscopic.
+    subroutine init_breit(j2max, specorbs)
         use parameter_def, only: NNNW
         use bcore_C, only: ICORE
         use bilst_C, only: FIRST, NTPI
         use decide_C, only: LTRANS
         use iounit_c, only: ISTDE
         use orb_C, only: NW, NCF
-        use wfac_C, only: WFACT
         use genintbreit1_I
         use genintbreit2_I
         use ichop_I
 
         integer, intent(in) :: j2max
+        logical, intent(in) :: specorbs(NNNW)
 
         ! These are the MPI parameters that need to be passed to different
         ! routines. We use the single core values.
         integer, parameter :: myid = 0, nprocs = 1
 
         integer :: i, j, N
 
+        ! Set up the the global array of specroscopic orbitals
+        breit_specorbs = specorbs
+
         ! We'll enable all parts of the Hamiltonian. E.g. AUXBLK relies on these
         ! flags to determine if certain things get initialized.
         LTRANS = .TRUE.
 
-        ! WFACT is the Breit scale factor. We set it to the default 1e-6
-        WFACT = 1d-6
-
         ! From rci3mpi_LINUX.f
         call genintbreit1(myid, nprocs, N, j2max)
         call genintbreit2(myid, nprocs, N, j2max)

src/matrixelements/rcisettings.f90

@@ -11,7 +11,8 @@ module grasp_rciqed_rcisettings
         ! Grid parameters
         !real(real64) :: RNT, H, N
         ! Hamiltonian parameters
-        logical :: breit_enabled, nms_enabled, sms_enabled, qed_vp_enabled
+        character(len=:), allocatable :: breit_mode
+        logical :: nms_enabled, sms_enabled, qed_vp_enabled
         ! For self-energy, we store:
         integer :: qed_se, qed_se_hydrogenic_cutoff
     end type rcisettings
@@ -59,7 +60,7 @@ function read_settings_toml(jobname, settings)
 
         ! TODO: Load nuclear model and GRID too!
 
-        call get_logical_default(table, "hamiltonian.breit", .false., settings%breit_enabled)
+        call get_string_default(table, "hamiltonian.breit", "n", settings%breit_mode)
         call get_logical_default(table, "hamiltonian.nms", .false., settings%nms_enabled)
         call get_logical_default(table, "hamiltonian.sms", .false., settings%sms_enabled)
         call get_logical_default(table, "hamiltonian.qed_vp", .false., settings%qed_vp_enabled)
@@ -89,7 +90,7 @@ end function read_settings_toml
 
     !> Writes a $(jobname).settings.toml file with the nuclear, grid and
     !! Hamiltonian settings that were used for the RCI run.
-    subroutine write_settings_toml(jobname, setype)
+    subroutine write_settings_toml(jobname, breit_mode, setype)
         use decide_C, only: LTRANS, LNMS, LSMS, LVP, LSE
         use def_C, only: Z, EMN, AUMAMU, FMTOAU
         use grid_C, only: RNT, H, N
@@ -98,10 +99,10 @@ subroutine write_settings_toml(jobname, setype)
         use grasp_rciqed_qed, only: nsetypes, setypes_short
 
         character(len=*), intent(in) :: jobname
-        ! TODO: passing setype explicitly here like this, while otherwise we
-        ! rely on the global common block modules is a bit of an inconsistent
-        ! hack to resolve a circular dependency problem. This should be fixed
-        ! at some point.
+        ! TODO: passing setype and breit_mopde explicitly here like this, even though
+        ! otherwise we rely on the global common block modules is a bit of an inconsistent
+        ! hack to resolve a circular dependency problem. This should be fixed at some point.
+        character(len=*), intent(in) :: breit_mode
         integer, intent(in) :: setype
 
         character(len=:), allocatable :: tomlfile
@@ -142,18 +143,26 @@ subroutine write_settings_toml(jobname, setype)
 
         write(toml_unit, '(a)') "[hamiltonian]"
         write(toml_unit, '(a)') "  # Contains the following booleans:"
-        write(toml_unit, '(a)') "  #   breit, nms, sms, qed_vp"
+        write(toml_unit, '(a)') "  #   qed_vp, nms, sms"
         write(toml_unit, '(a)') "  # Each indicates if the corresponding part of the Hamiltonian"
         write(toml_unit, '(a)') "  # was enabled. If it is omitted, it is assumed to have been off."
         write(toml_unit, '(a)') "  #"
+        write(toml_unit, '(a)') "  # Can also contain breit (string) which indicates the mode of the"
+        write(toml_unit, '(a)') "  # transverse photon operator. Possible values are 'breit', 'specorbs'"
+        write(toml_unit, '(a)') "  # and 'full'. If set to 'n', the operator is completely disabled."
+        write(toml_unit, '(a)') "  #"
         write(toml_unit, '(a)') "  # May also contain qed_se (string), which indicates that a particular"
         write(toml_unit, '(a)') "  # QED self-energy operator was also included in the Hamiltonian."
         write(toml_unit, '(a)') "  # The possible values are: 'hydrogenic', 'qedmod', 'flambaum' and 'pyykkoe'"
         write(toml_unit, '(a)') "  #"
         write(toml_unit, '(a)') "  # Finally, for the hydrogenic QED self-energy, qed_se_hydrogenic_cutoff"
         write(toml_unit, '(a)') "  # (integer) may be defined, which sets the n-quantum number cutoff. If not"
         write(toml_unit, '(a)') "  # present, the implementation default is used."
-        if(LTRANS) write(toml_unit, '(a)') "  breit = true"
+        if(LTRANS) then
+            write(toml_unit, '("  breit = """,a,"""")') trim(breit_mode)
+        else
+            write(toml_unit, '("  breit = ""n""")')
+        endif
         if(LNMS) write(toml_unit, '(a)') "  nms = true"
         if(LSMS) write(toml_unit, '(a)') "  sms = true"
         if(LVP) write(toml_unit, '(a)') "  qed_vp = true"

src/pt.f90

@@ -8,6 +8,7 @@
 !!
 program rci_qed_pt
     use, intrinsic :: iso_fortran_env, only: real64, dp => real64
+    use parameter_def, only: NNNW
     use grasp_lib9290, only: init_isocw_full
     use grasp_rciqed_lib9290_files, only: load_mixing
     use grasp_rciqed_lib9290_csls, only: ncsfs_global
@@ -120,7 +121,9 @@ program rci_qed_pt
     call load_mixing(file_mixing)
     print '(a)', ">>> INITIALIZING rci commons"
     call init_rkintc(j2max)
-    call init_breit(j2max)
+    ! TODO: FIX THIS -- should be generic
+    print '(a)', ">>> INITIALIZING fully frequency-dependent Breit"
+    call init_breit(j2max, (/ (.true., k=1,NNNW) /))
     call qedvp_init
     call init_mass_shifts
 
@@ -141,7 +144,11 @@ program rci_qed_pt
 
     print "(a)", "Hamiltonian parts enabled in RCI calculation:"
     print '(" ",a,"  ",a)', check(.true.), "Dirac + nuclear potential"
-    print '(" ",a,"  ",a)', check(settings%breit_enabled), "Breit"
+    if(settings%breit_mode == "n") then
+        print '(" ",a,"  ",a)', check(.false.), "Breit"
+    else
+        print '(" ",a,"  Breit (mode: ",a,")")', check(.true.), settings%breit_mode
+    endif
     print '(" ",a,"  ",a)', check(settings%nms_enabled), "Normal mass shift"
     print '(" ",a,"  ",a)', check(settings%sms_enabled), "Special mass shift"
     print '(" ",a,"  ",a)', check(settings%qed_vp_enabled), "QED vacuum polarization"
@@ -158,11 +165,12 @@ program rci_qed_pt
 
     hcs_cat(1) = .true.
     hcs_cat(2) = .true.
-    hcs_cat(3) = settings%breit_enabled
+    ! TODO: actually, we should split up by contribution?
+    hcs_cat(3) = (settings%breit_mode /= "n")
     hcs_cat(4) = settings%nms_enabled
     hcs_cat(5) = settings%sms_enabled
     hcs_cat(6) = settings%qed_vp_enabled
-    hcs_cat(7) = .not.(settings%qed_se == 0)
+    hcs_cat(7) = (settings%qed_se /= 0)
 
     ! Also set the hydrogenic settings in qedcut_C appropriately. LSE needs to be set to
     ! .true. since QED_SLFEN checks for it when applying NQEDMAX for some reason.

src/rci-qed.f90

@@ -68,7 +68,7 @@ PROGRAM RCI90mpi
       USE strsum_I
       USE factt_I
       USE matrix_I
-      use grasp_rciqed, only: IMCDF => res_unit, setype
+      use grasp_rciqed, only: IMCDF => res_unit, setype, breit_mode
       use grasp_rciqed_rcisettings, only: write_settings_toml
       IMPLICIT NONE
 !-----------------------------------------------
@@ -219,7 +219,7 @@ PROGRAM RCI90mpi
          GO TO 99
       ENDIF
 
-      if(myid == 0) call write_settings_toml(name, setype)
+      if(myid == 0) call write_settings_toml(name, breit_mode, setype)
 !=======================================================================
 !  Execution finished; Statistics output
 !=======================================================================

src/rci/getcid.f90

@@ -32,12 +32,12 @@ SUBROUTINE GETCID (isofile, rwffile, idblk)
       USE nsmdat_C
       USE orb_C
       USE wave_C
-      USE wfac_C
       USE blim_C
       USE qedcut_C
       USE mpi_C
-      use getcid_I, only: getcid_qedse
-      use grasp_rciqed, only: IMCDF => res_unit, setype
+      use getcid_I, only: getcid_specorbs, getcid_breit, getcid_qedse
+      use grasp_rciqed, only: IMCDF => res_unit, setype, isspecorb
+      use grasp_rciqed_breit, only: breit_specorbs
 !-----------------------------------------------
 !   I n t e r f a c e   B l o c k s
 !-----------------------------------------------
@@ -138,6 +138,10 @@ SUBROUTINE GETCID (isofile, rwffile, idblk)
                                                MPI_COMM_WORLD,ierr)
       ENDIF ! .NOT. LFORDR
 
+      ! Determine spectroscopic orbitals
+      if (myid == 0) call getcid_specorbs
+      CALL MPI_Bcast (isspecorb, size(isspecorb), MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+
 !*****************************************************************
 !
 ! Pre-run ?
@@ -166,24 +170,12 @@ SUBROUTINE GETCID (isofile, rwffile, idblk)
 !     ENDIF
 !*****************************************************************
 !
-! Include transverse ?
-! Quantities to be obtained: LTRANS, WFACT
-!
-      IF (myid .EQ. 0) THEN
-         WRITE (istde,*) 'Include contribution of H (Transverse)?'
-         LTRANS = GETYN ()
-         WRITE (istde,*) 'Modify all transverse photon frequencies?'
-         YES = GETYN ()
-         IF (YES) THEN
-            WRITE (istde,*) 'Enter the scale factor:'
-            READ *, WFACT
-         ELSE
-            WFACT = 1.0D00
-         ENDIF
-      ENDIF ! myid .EQ. 0
+
+      ! Include Breit? Also broadcast the settings to other nodes.
+      if (myid == 0) call getcid_breit
       CALL MPI_Bcast (LTRANS, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
-      CALL MPI_Bcast (WFACT, 1, MPI_DOUBLE_PRECISION, 0,            &
-                           MPI_COMM_WORLD, ierr)
+      CALL MPI_Bcast (breit_specorbs, size(breit_specorbs), MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+
 !
 ! Other interactions ? One logical for each case. Done altogether
 !
@@ -204,12 +196,12 @@ SUBROUTINE GETCID (isofile, rwffile, idblk)
          ELSE
           WRITE(734,'(a)')'n            ! Contribution of H Transverse?'
          END IF
-         IF (YES) THEN
-           WRITE(734,'(a)') 'y            ! Modify photon frequencies?'
-           WRITE(734,*) WFACT,'! Scale factor'
-         ELSE
-           WRITE(734,'(a)') 'n            ! Modify photon frequencies?'
-         END IF
+         !IF (YES) THEN
+         !  WRITE(734,'(a)') 'y            ! Modify photon frequencies?'
+         !  WRITE(734,*) WFACT,'! Scale factor'
+         !ELSE
+         !  WRITE(734,'(a)') 'n            ! Modify photon frequencies?'
+         !END IF
          IF (LVP) THEN
            WRITE(734,'(a)') 'y            ! Vacuum polarization?'
          ELSE
@@ -332,7 +324,8 @@ SUBROUTINE GETCID (isofile, rwffile, idblk)
 !   This is the sixth record on the file
 !
       WRITE (imcdf) C, LFORDR, (ICCUTblk(i), i = 1, nblock),           &
-                    LTRANS, WFACT, LVP, LNMS, LSMS
+                    LTRANS, LVP, LNMS, LSMS
+      WRITE (imcdf) (breit_specorbs(i), i = 1, NW)
 !
 !   Write the grid data to the  .res  file; this is the seventh
 !   record on the file