module nmr_shielding_mod

  implicit none

  character(len=*), parameter :: module_name = "nmr_shielding_mod"

  private
  public nmr_shielding

contains

  subroutine nmr_shielding_C(c_handle) bind(C, name="nmr_shielding")
    use c_interop, only: oqp_handle_t, oqp_handle_get_info
    use types, only: information
    type(oqp_handle_t) :: c_handle
    type(information), pointer :: inf
    inf => oqp_handle_get_info(c_handle)
    call nmr_shielding(inf)
  end subroutine nmr_shielding_C

!> @brief NMR nuclear magnetic shielding tensors.
!> @details Current validated scope: RHF and closed-shell DFT, common gauge
!>  origin (CGO). The user-facing Python layer recognizes `properties.nmr_gauge`
!>  with CGO as the default and GIAO as a gated development option; GIAO NMR
!>  shielding is not yet validated and does not enter this Fortran CGO pathway.
!>  the uncoupled and the coupled (CPHF/CPKS) paramagnetic responses are reported.
!>  The coupled response includes the exact-exchange response of the imaginary
!>  antisymmetric first-order magnetic density (Phase 0); the Coulomb and
!>  semi-local XC-kernel responses vanish by symmetry, so the coupling is exact
!>  exchange scaled by the exchange fraction c_x (zero for pure functionals, where
!>  coupled == uncoupled). The isotropic shielding is sigma = sigma_dia +
!>  sigma_para per nucleus.
!>
!> Phase-0 validation (H2O/STO-3G, CGO at COM; common-gauge reference values in
!> tests/fixtures/nmr/cgo_reference.json):
!>   - HF coupled para matches the oracle exactly (O -230.63, H 3.506 ppm).
!>   - Coulomb response of P^B ~0; exact-exchange response nonzero (gates 1-2).
!>   - Pure PBE: coupled == uncoupled (gate 3). HF/hybrid coupled != uncoupled,
!>     with the coupling scaling with c_x (gates 4, 6).
!>
!> Validation (H2O/STO-3G, RHF, CGO at the center of mass; reference values):
!>   - Diamagnetic term matches the standard common-gauge diamagnetic to ~6 significant
!>     figures (O 411.418, H 28.062 ppm).
!>   - Paramagnetic term (MO transform of the orbital-Zeeman and PSO operators,
!>     occupied-virtual sum-over-states, 2*alpha^2 prefactor) matches the reference
!>     uncoupled reference for BOTH atoms (O para -113.63, H para 1.785 ppm;
!>     totals 297.79 / 29.85 ppm).
!>   - The PSO operator is anti-Hermitian; `pso_integrals` returns it exactly
!>     antisymmetric (max|diag| and max|A+A^T| are reported below as a check).
  subroutine nmr_shielding(infos)
    use io_constants, only: iw
    use oqp_tagarray_driver
    use basis_tools, only: basis_set
    use messages, only: show_message, with_abort
    use types, only: information
    use int1, only: angular_momentum_integrals, nmr_dia_shielding, pso_integrals

    implicit none

    character(len=*), parameter :: subroutine_name = "nmr_shielding"
    ! CODATA fine-structure constant and derived prefactors
    real(kind=8), parameter :: ALPHA = 1.0d0/137.035999084d0
    real(kind=8), parameter :: halfa2 = 0.5d0*ALPHA*ALPHA
    real(kind=8), parameter :: twoa2 = 2.0d0*ALPHA*ALPHA
    real(kind=8), parameter :: PPM = 1.0d6

    type(information), target, intent(inout) :: infos

    integer :: nbf, nbf2, ok
    logical :: urohf
    type(basis_set), pointer :: basis
    real(kind=8), allocatable :: amom(:,:)      ! packed Lx,Ly,Lz (lower triangle)
    real(kind=8), allocatable :: lfull(:,:,:)   ! full antisymmetric (nbf,nbf,3)
    real(kind=8), allocatable :: gdia(:,:,:)    ! diamagnetic contracted integrals (3,3,nat)
    real(kind=8), allocatable :: sig_dia(:,:,:) ! diamagnetic shielding tensor (3,3,nat)
    real(kind=8), allocatable :: coords(:,:)    ! nuclear coordinates (3,nat)
    real(kind=8), allocatable :: siso_dia(:)    ! isotropic diamagnetic shielding (ppm)
    real(kind=8), allocatable :: pso_full(:,:,:)! full antisymmetric PSO (nbf,nbf,3)
    real(kind=8), allocatable :: moL(:,:,:)     ! orbital Zeeman in MO basis (nmo,nmo,3)
    real(kind=8), allocatable :: moP(:,:,:)     ! PSO in MO basis (nmo,nmo,3)
    real(kind=8), allocatable :: sig_para(:,:,:)! paramagnetic shielding tensor (3,3,nat)
    real(kind=8), allocatable :: siso_para(:), siso_tot(:)
    ! Phase 0: coupled (CPHF/CPKS) magnetic response
    real(kind=8), allocatable :: rcoup(:,:)     ! coupled response vectors (lvir,3)
    real(kind=8), allocatable :: sig_para_c(:,:,:), siso_para_c(:), siso_tot_c(:)
    real(kind=8) :: scale_exch, pb_asym, jnorm, knorm
    integer :: nvir, lvir
    logical :: is_dft
    real(kind=8) :: o(3), com(3), trg, pso_diag_max, pso_asym_max
    integer :: nat, i, c, t, s, nocc, nmo

    real(kind=8), contiguous, pointer :: dmat_a(:)
    real(kind=8), contiguous, pointer :: mo_a(:,:)
    real(kind=8), contiguous, pointer :: mo_e_a(:)
    integer(4) :: status

    urohf = infos%control%scftype == 2 .or. infos%control%scftype == 3

    open (unit=IW, file=infos%log_filename, position="append")

    basis => infos%basis
    basis%atoms => infos%atoms

    nbf = basis%nbf
    nbf2 = nbf*(nbf+1)/2

    write(iw,'(2/)')
    write(iw,'(4x,a)') '======================================'
    write(iw,'(4x,a)') 'NMR nuclear magnetic shielding (CGO)'
    write(iw,'(4x,a)') '======================================'
    write(iw,'(4x,a)') 'Gauge formulation: CGO (common gauge origin). For gauge-origin-'// &
      'independent results use properties.nmr_gauge=giao.'
    call flush(iw)

    if (urohf) then
      call show_message('NMR shielding currently supports closed-shell &
        &(RHF / pure-DFT) references only', with_abort)
    end if

    ! Not-implemented classes must abort instead of silently producing wrong
    ! shieldings:
    !  - CAM/range-separated hybrids: the coupled magnetic response and the
    !    two-electron derivative use the global exchange fraction only; the
    !    range-separation attenuation (alpha/beta/mu) is not wired in.
    !  - meta-GGAs: the tau channel of the XC response is not implemented.
    !  - ECP: no effective-core magnetic-derivative term is implemented.
    if (infos%control%hamilton == 20) then
      if (infos%dft%cam_flag) then
        call show_message('NMR shielding with range-separated (CAM) functionals &
          &is not implemented', with_abort)
      end if
      if (infos%functional%needtau) then
        call show_message('NMR shielding with meta-GGA (tau-dependent) functionals &
          &is not implemented', with_abort)
      end if
    end if
    if (allocated(infos%basis%ecp_zn_num)) then
      if (any(infos%basis%ecp_zn_num /= 0)) then
        call show_message('NMR shielding with ECP basis sets is not implemented', &
          with_abort)
      end if
    end if

    ! Confirm the SCF density is present (used by later stages)
    call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a, status)
    call check_status(status, module_name, subroutine_name, OQP_DM_A)

    ! Gauge origin: center of mass (default)
    nat = ubound(basis%atoms%zn,1)
    com = 0
    do i = 1, nat
      com = com + basis%atoms%xyz(:,i)*basis%atoms%mass(i)
    end do
    com = com / sum(basis%atoms%mass)
    o = com

    write(iw,'(/4x,a)') 'Gauge origin (Bohr):'
    write(iw,'(4x,a,3f15.8)') 'O = ', o

    ! Angular momentum integrals about the gauge origin (packed, antisymmetric)
    allocate(amom(nbf2,3), source=0.0d0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)

    call angular_momentum_integrals(basis, amom, o)

    ! Expand each component to a full antisymmetric nbf x nbf matrix (the
    ! orbital-Zeeman / magnetic-field perturbation used by the paramagnetic term).
    allocate(lfull(nbf,nbf,3), source=0.0d0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    do c = 1, 3
      call expand_antisym(amom(:,c), lfull(:,:,c), nbf)
    end do

    !-------------------------------------------------------------------------
    ! Diamagnetic term
    !   sigma^dia_{ts}(N) = (alpha^2/2) [ delta_ts*Tr(g^N) - g^N_{s,t} ]
    ! where g^N_{ab} = sum_{mu,nu} P_{mu,nu} <mu|(r-O)_a (r-R_N)_b/|r-R_N|^3|nu>
    !-------------------------------------------------------------------------
    allocate(gdia(3,3,nat), coords(3,nat), source=0.0d0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    do i = 1, nat
      coords(:,i) = basis%atoms%xyz(:,i)
    end do

    call nmr_dia_shielding(basis, dmat_a, o, coords, nat, gdia)

    allocate(sig_dia(3,3,nat), siso_dia(nat), source=0.0d0)
    do i = 1, nat
      trg = gdia(1,1,i) + gdia(2,2,i) + gdia(3,3,i)
      do t = 1, 3
        do s = 1, 3
          sig_dia(t,s,i) = halfa2 * (merge(trg, 0.0d0, t==s) - gdia(s,t,i))
        end do
      end do
      siso_dia(i) = (sig_dia(1,1,i)+sig_dia(2,2,i)+sig_dia(3,3,i))/3.0d0 * PPM
    end do

    !-------------------------------------------------------------------------
    ! Paramagnetic term (uncoupled; exact CPKS for pure functionals/HF-uncoupled)
    !   sigma^para_{xy}(N) = 2 alpha^2 sum_{i occ, a vir}
    !                          MO_x(a,i) * PSO^N_y(i,a) / (eps_a - eps_i)
    ! MO_x = C^T A_O[x] C  (orbital Zeeman, angular momentum about O)
    ! PSO^N_y = C^T A_PSO^N[y] C
    !-------------------------------------------------------------------------
    call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a, status)
    call check_status(status, module_name, subroutine_name, OQP_VEC_MO_A)
    call tagarray_get_data(infos%dat, OQP_E_MO_A, mo_e_a, status)
    call check_status(status, module_name, subroutine_name, OQP_E_MO_A)

    nocc = int(infos%mol_prop%nocc)
    nmo  = size(mo_e_a)

    ! Orbital-Zeeman operator in MO basis (3 components)
    allocate(moL(nmo,nmo,3), source=0.0d0)
    do c = 1, 3
      call ao_to_mo(lfull(:,:,c), mo_a, moL(:,:,c), nbf, nmo)
    end do

    ! ----------------------------------------------------------------------
    ! Phase 0: coupled magnetic response (CPHF/CPKS).
    ! The first-order magnetic density is imaginary/antisymmetric, so the
    ! Coulomb and (semi-local) XC-kernel responses vanish; only the exact
    ! exchange response survives, scaled by the exchange fraction c_x. For
    ! pure functionals (c_x = 0) the coupled response equals the uncoupled one.
    ! ----------------------------------------------------------------------
    nvir = nmo - nocc
    lvir = nocc*nvir
    is_dft = infos%control%hamilton == 20
    scale_exch = 1.0d0
    if (is_dft) scale_exch = infos%dft%HFscale
    allocate(rcoup(lvir,3), source=0.0d0)
    call compute_coupled_para(infos, basis, mo_a, mo_e_a, nocc, nbf, nvir, &
                              moL, scale_exch, rcoup, pb_asym, jnorm, knorm)

    allocate(pso_full(nbf,nbf,3), moP(nmo,nmo,3), source=0.0d0)
    allocate(sig_para(3,3,nat), siso_para(nat), siso_tot(nat), source=0.0d0)
    allocate(sig_para_c(3,3,nat), siso_para_c(nat), siso_tot_c(nat), source=0.0d0)

    pso_diag_max = 0.0d0
    pso_asym_max = 0.0d0
    do i = 1, nat
      ! Full antisymmetric PSO matrices A_a = [(r-R_N) x grad]_a/|r-R_N|^3
      call pso_integrals(basis, coords(:,i), pso_full)
      do c = 1, 3
        ! Diagnostics: the PSO operator is anti-Hermitian, so the diagonal and
        ! the symmetric part must vanish.
        do t = 1, nbf
          pso_diag_max = max(pso_diag_max, abs(pso_full(t,t,c)))
          do s = 1, nbf
            pso_asym_max = max(pso_asym_max, abs(pso_full(t,s,c)+pso_full(s,t,c)))
          end do
        end do
        call ao_to_mo(pso_full(:,:,c), mo_a, moP(:,:,c), nbf, nmo)
      end do
      do t = 1, 3
        do s = 1, 3
          sig_para(t,s,i)   = twoa2 * sum_ov(moL(:,:,t), moP(:,:,s), mo_e_a, nocc, nmo)
          sig_para_c(t,s,i) = -twoa2 * sum_resp(rcoup(:,t), moP(:,:,s), nocc, nvir)
        end do
      end do
      siso_para(i)   = (sig_para(1,1,i)+sig_para(2,2,i)+sig_para(3,3,i))/3.0d0 * PPM
      siso_tot(i)    = siso_dia(i) + siso_para(i)
      siso_para_c(i) = (sig_para_c(1,1,i)+sig_para_c(2,2,i)+sig_para_c(3,3,i))/3.0d0 * PPM
      siso_tot_c(i)  = siso_dia(i) + siso_para_c(i)
    end do

    write(iw,'(/4x,a,f8.4,a)') 'Isotropic shielding (CGO, ppm)   [exact-exchange c_x =', &
           scale_exch, ']'
    write(iw,'(4x,a)')  '   Atom    Z   sigma_dia   para_uncoupled   para_coupled'// &
           '   total_uncoupled   total_coupled'
    do i = 1, nat
      write(iw,'(4x,i6,f6.1,5f16.6)') i, basis%atoms%zn(i), &
             siso_dia(i), siso_para(i), siso_para_c(i), siso_tot(i), siso_tot_c(i)
    end do

    ! ---- Store isotropic shielding (ppm) to a tagarray for JSON output ----
    block
      real(kind=8), contiguous, pointer :: nmrout(:)
      call infos%dat%alloc_or_die(OQP_nmr_shielding, (/ 5*nat /), nmrout, &
                                  description=OQP_nmr_shielding_comment)
      do i = 1, nat
        nmrout(5*(i-1)+1) = siso_dia(i)
        nmrout(5*(i-1)+2) = siso_para(i)
        nmrout(5*(i-1)+3) = siso_para_c(i)
        nmrout(5*(i-1)+4) = siso_tot(i)
        nmrout(5*(i-1)+5) = siso_tot_c(i)
      end do
    end block

    ! ---- Phase-0 validation gates (reported as diagnostics) ----
    write(iw,'(/4x,a)') 'Phase-0 magnetic-response gates:'
    write(iw,'(4x,a,es12.3)') '  gate0  max|P^B + (P^B)^T|        = ', pb_asym
    write(iw,'(4x,a,es12.3)') '  gate1  ||J(P^B)|| (Coulomb)      = ', jnorm
    write(iw,'(4x,a,es12.3)') '  gate2  ||K(P^B)|| (exact exch.)  = ', knorm
    write(iw,'(4x,a,2es12.3)') '  PSO    max|diag|, max|A+A^T|    = ', &
           pso_diag_max, pso_asym_max
    call flush(iw)

    deallocate(amom, lfull, gdia, coords, sig_dia, siso_dia)
    deallocate(pso_full, moL, moP, sig_para, siso_para, siso_tot)
    deallocate(rcoup, sig_para_c, siso_para_c, siso_tot_c)
    close(iw)

  end subroutine nmr_shielding

!> @brief Expand a packed lower-triangular antisymmetric matrix to full form.
!> @details Packed storage holds the bra>=ket elements A(p,q) (p>=q) at index
!>  q + p*(p-1)/2. The full matrix satisfies A(q,p) = -A(p,q), zero diagonal.
  subroutine expand_antisym(packed, full, n)
    real(kind=8), intent(in)  :: packed(:)
    real(kind=8), intent(out) :: full(:,:)
    integer, intent(in) :: n
    integer :: p, q, idx
    full = 0.0d0
    do p = 1, n
      do q = 1, p
        idx = q + p*(p-1)/2
        full(p,q) =  packed(idx)
        full(q,p) = -packed(idx)
      end do
    end do
  end subroutine expand_antisym

!> @brief Transform a full AO matrix to the MO basis: M = C^T A C.
  subroutine ao_to_mo(a_ao, c, m_mo, nbf, nmo)
    real(kind=8), intent(in)  :: a_ao(:,:)   ! (nbf,nbf)
    real(kind=8), intent(in)  :: c(:,:)      ! (nbf,nmo)
    real(kind=8), intent(out) :: m_mo(:,:)   ! (nmo,nmo)
    integer, intent(in) :: nbf, nmo
    real(kind=8), allocatable :: tmp(:,:)
    allocate(tmp(nbf,nmo))
    tmp = matmul(a_ao, c(:,1:nmo))
    m_mo = matmul(transpose(c(:,1:nmo)), tmp)
    deallocate(tmp)
  end subroutine ao_to_mo

!> @brief Occupied-virtual sum-over-states contraction
!>   sum_{i occ, a vir} L(a,i) * P(i,a) / (eps_a - eps_i)
  function sum_ov(lmo, pmo, e, nocc, nmo) result(val)
    real(kind=8), intent(in) :: lmo(:,:), pmo(:,:), e(:)
    integer, intent(in) :: nocc, nmo
    real(kind=8) :: val
    integer :: i, a
    val = 0.0d0
    do i = 1, nocc
      do a = nocc+1, nmo
        val = val + lmo(a,i)*pmo(i,a)/(e(a)-e(i))
      end do
    end do
  end function sum_ov

!> @brief Contract a coupled response vector (occ-vir, length nocc*nvir, packed
!>  as k=(a-1)*nocc+i to match iatogen) with the PSO occ-vir block.
  function sum_resp(rt, pmo, nocc, nvir) result(val)
    real(kind=8), intent(in) :: rt(:)        ! (nocc*nvir)
    real(kind=8), intent(in) :: pmo(:,:)     ! PSO[s] in MO basis (nmo,nmo)
    integer, intent(in) :: nocc, nvir
    real(kind=8) :: val
    integer :: i, a, k
    val = 0.0d0
    k = 0
    do a = 1, nvir
      do i = 1, nocc
        k = k + 1
        val = val + rt(k)*pmo(i, nocc+a)
      end do
    end do
  end function sum_resp

!> @brief Build the antisymmetric AO first-order magnetic density from an
!>  occ-vir response vector: pa = C * (av - av^T) * C^T, av(occ,vir) = rin.
  subroutine magnetic_pb_density(mo, pa, av, nbf, nocc, rin)
    use tdhf_lib, only: iatogen
    use mathlib, only: orthogonal_transform
    real(kind=8), intent(in) :: mo(:,:)
    real(kind=8), intent(inout), target :: pa(:,:,:)
    real(kind=8), intent(inout) :: av(:,:)
    integer, intent(in) :: nbf, nocc
    real(kind=8), intent(in) :: rin(:)
    call iatogen(rin, av, nocc, nocc)
    av = av - transpose(av)
    call orthogonal_transform('t', nbf, mo, av, pa(:,:,1))
  end subroutine magnetic_pb_density

!> @brief Solve the coupled (CPHF/CPKS) magnetic response for the three field
!>  components and report the Phase-0 gate diagnostics.
!> @details Fixed-point solve of  (eps_a-eps_i) R + c_x*K[P^B(R)] = b,
!>  with b = orbital-Zeeman occ-vir block and K the exact-exchange image of the
!>  antisymmetric first-order density. Coulomb (J) and the semi-local XC kernel
!>  do not contribute (imaginary antisymmetric density). For c_x = 0 the loop is
!>  skipped and R = b/(eps_a-eps_i) (uncoupled).
  subroutine compute_coupled_para(infos, basis, mo, mo_e, nocc, nbf, nvir, &
                                  moL, scale_exch, rcoup, pb_asym, jnorm, knorm)
    use int2_compute, only: int2_compute_t
    use tdhf_lib, only: int2_td_data_t, mntoia
    use types, only: information
    use basis_tools, only: basis_set
    use messages, only: show_message
    implicit none
    type(information), target, intent(inout) :: infos
    type(basis_set), intent(in) :: basis
    real(kind=8), intent(in) :: mo(:,:), mo_e(:), moL(:,:,:)
    integer, intent(in) :: nocc, nbf, nvir
    real(kind=8), intent(in) :: scale_exch
    real(kind=8), intent(out) :: rcoup(:,:)         ! (lvir,3)
    real(kind=8), intent(out) :: pb_asym, jnorm, knorm

    integer, parameter :: maxit = 100
    real(kind=8), parameter :: tol = 1.0d-9, half = 0.5d0
    real(kind=8), parameter :: kappa = 1.0d0        ! coupling sign (validated vs the reference)
    type(int2_compute_t) :: int2_driver
    type(int2_td_data_t), target :: kdat, kdat1, jdat
    real(kind=8), allocatable, target :: pa(:,:,:)
    real(kind=8), allocatable :: av(:,:), bb(:,:), dd(:), gx(:), rprev(:), gao(:,:)
    integer :: lvir, t, i, a, k, it
    logical :: conv

    lvir = nocc*nvir
    allocate(pa(nbf,nbf,1), av(nbf,nbf), gao(nbf,nbf), &
             bb(lvir,3), dd(lvir), gx(lvir), rprev(lvir), source=0.0d0)

    ! RHS (orbital-Zeeman occ-vir block) and orbital-energy denominators
    do t = 1, 3
      k = 0
      do a = 1, nvir
        do i = 1, nocc
          k = k + 1
          bb(k,t) = moL(i, nocc+a, t)
        end do
      end do
    end do
    k = 0
    do a = 1, nvir
      do i = 1, nocc
        k = k + 1
        dd(k) = mo_e(nocc+a) - mo_e(i)
      end do
    end do

    do t = 1, 3
      rcoup(:,t) = bb(:,t)/dd                        ! uncoupled start
    end do

    call int2_driver%init(basis, infos)
    ! NMR uses the native Rys ERI path only; never route through libint.
    int2_driver%rys_only = .true.
    call int2_driver%set_screening()

    ! Coupled iteration (skipped for pure functionals, c_x = 0)
    if (abs(scale_exch) > 1.0d-12) then
      kdat = int2_td_data_t(d2=pa, int_apb=.false., int_amb=.true., &
                            tamm_dancoff=.false., scale_exchange=scale_exch)
      do t = 1, 3
        conv = .false.
        do it = 1, maxit
          rprev = rcoup(:,t)
          call magnetic_pb_density(mo, pa, av, nbf, nocc, rcoup(:,t))
          call int2_driver%run(kdat)
          gao = half*kdat%amb(:,:,1,1)
          call mntoia(gao, gx, mo, mo, nocc, nocc)
          rcoup(:,t) = (bb(:,t) - kappa*gx)/dd
          if (maxval(abs(rcoup(:,t)-rprev)) < tol) then
            conv = .true.
            exit
          end if
        end do
        if (.not. conv) then
          call show_message('WARNING: NMR coupled magnetic response (CPHF) did &
            &not converge within the iteration limit; shieldings may be inaccurate')
        end if
      end do
    end if

    ! ---- Gate diagnostics on the converged z-component first-order density ----
    call magnetic_pb_density(mo, pa, av, nbf, nocc, rcoup(:,3))
    pb_asym = maxval(abs(pa(:,:,1) + transpose(pa(:,:,1))))      ! gate 0

    jdat = int2_td_data_t(d2=pa, int_apb=.true., int_amb=.false., &
                          tamm_dancoff=.false., scale_exchange=0.0d0)
    call int2_driver%run(jdat)
    jnorm = sqrt(sum((half*jdat%apb(:,:,1,1))**2))               ! gate 1 (Coulomb)

    kdat1 = int2_td_data_t(d2=pa, int_apb=.false., int_amb=.true., &
                           tamm_dancoff=.false., scale_exchange=1.0d0)
    call int2_driver%run(kdat1)
    knorm = sqrt(sum((half*kdat1%amb(:,:,1,1))**2))              ! gate 2 (exchange)

    deallocate(pa, av, gao, bb, dd, gx, rprev)
  end subroutine compute_coupled_para

end module nmr_shielding_mod
