#include "../include/assert.fh"
module xray_fourier_module
! This module has the basic routines for calculating X-ray forces.
! Fourier (reciprocal) space values are stored as list of H,K,L values rather
! than a 3D array. Normally, only H,K,L indices with an observed value are
! saved.
!
! These routines all pass F90-style dimensioned arrays.
! This version has several dependent modules commented out using CPP
! conditionals.
!
!  SUBROUTINES:
!
!  fourier_dXYZBQ_dF  --  Calculate the derivative of the coordinate,
!                              B, or occupancy versus the structure factor.
!                              This one uses the direct (exact) Fourier.
!
!  FFT_dXYZBQ_dF      --  Same as above, using the FFT.
!
!  dF_dTarget         --  Calculate a structure-factor restraint force
!                              from the difference of Fobs and Fcalc. The
!                              residual restraint target-function is essentially a
!                              harmonic restraint in reciprocal space.
!
!  FFT_Fcalc          --  Calculate structure factors from coordinates.
!
!  CNS_FFTab               --  This reduces the FFT grid to the observerved HKL list.
!                              and therefore defines the choice of reciprocal asymmetric
!                              unit. It was derived from CNS. This is actually
!                              somewhat complicated when symmetry is present, due to
!                              many special-case tests for reflections tha lie
!                              on symmetry planes or axes.
!
! FUNCTIONS:
!
! atom_scatter_factor_mss4 --  Calculate the atomic scatter (f) at a specific resolution, given
!                              given the Gaussian coefficients, and a modified resolution
!                              defined as -S*S/4.0, where S is the reciprocal
!                              resolution. (mss4 stands for "Minus S Squared over 4")
!                              See comments at the beginning of xray_fourier_Fcalc()
   use xray_globals_module
   use constants, only: M_TWOPI => TWOPI
   use xray_reciprocal_space_module
   use xray_real_space_module
   use xray_fftpack_module
   implicit none

   !-------------------------------------------------------------------
   ! Descriptor for the FFT implementation:
   type(fft3c_descriptor), save :: fft3_desc

   !-------------------------------------------------------------------
contains

   !-------------------------------------------------------------------
   ! Thus is a generic front-end to utilize a specific FFT library.
   subroutine FFT_setup()
      integer :: i
      if (spacegroup_number/=1) stop 'ERROR: only P1 symmetry is supported.'
      write(stdout,'(A)') '---------------- XRAY FFT SETUP -------------------'
      write(stdout,'(A,3F9.3,3F7.2)') 'Unit cell: ',unit_cell
      write(stdout,'(A,F9.5)') 'FFT: Requested grid spacing: ',fft_grid_spacing
      grid_size=ceiling(unit_cell(1:3)/fft_grid_spacing)
      write(stdout,'(A,3I8)') 'FFT: Initial grid size: ',grid_size
      do i=1,3
         call fft_adjust_size(grid_size(i),2)
      end do
      write(stdout,'(A,3I8)') 'FFT: Adjusted grid size: ',grid_size
      frac_to_grid = grid_size
      write(stdout,*) frac_to_grid
      angstrom_to_grid = frac_to_grid / unit_cell(1:3)
      write(stdout,*) angstrom_to_grid
      do i=1,3
        orth_to_grid(:,i) = orth_to_frac(:,i) * frac_to_grid(i)
        grid_to_orth(:,i) = frac_to_orth(:,i) / frac_to_grid(i)
      end do
      write(stdout,'(A,3F8.3)') 'FFT: Actual grid spacing: ',angstrom_to_grid**(-1)
      call fft3c_new(fft3_desc,grid_size)
      allocate(density_map(0:grid_size(1)-1, &
                           0:grid_size(2)-1, &
                           0:grid_size(3)))
   end subroutine FFT_setup

   subroutine FFT_forward()
   end subroutine FFT_forward

   subroutine FFT_backward()
   end subroutine FFT_backward

   !-------------------------------------------------------------------
   ! Caution: Some literature uses S to represent S^2
   !
   ! S == d* = sqrt(sum(HKL * orth_to_frac)^2) = sqrt(-4*mSS4)
   ! mSS4 = -S*S/4
   ! mSS4 is stored instead of S because it is more relevant to the formulas used.

   subroutine fourier_Fcalc( &
            num_hkl,hkl,Fcalc,mSS4,hkl_selected, & ! reflections
      num_atoms,xyz,tempFactor,occupancy,scatter_type_index)
      ! Globals: scatter_cofficients
      implicit none
      integer, intent(in) :: num_hkl
      integer, intent(in) :: hkl(3,num_hkl)
      complex(real_kind), intent(out) :: Fcalc(num_hkl)       ! Fcalc (from the atomic coordinates)
      real(real_kind), intent(in) :: mSS4(num_hkl)
      integer, intent(in), optional :: hkl_selected(num_hkl)  ! Subset of reflections to use in this calculation (default = all)

      integer, intent(in) :: num_atoms
      real(real_kind), intent(in), target :: xyz(3,num_atoms)               ! (3,n) atomic coordinates
      real(real_kind), intent(in) :: tempFactor(num_atoms)                  ! atomic temperature (B) factors
      real(real_kind), intent(in), target, optional :: occupancy(num_atoms) ! atomic occupancy (default = 1.0)
      integer, intent(in), target :: scatter_type_index(num_atoms) ! index into the scatter - type coeffs table for each atom
      ! locals
      ! integer :: work_num_atoms
      integer :: ihkl, i, alloc_status
      ! automatic
      real(real_kind) :: atomic_scatter_factor(num_scatter_types)
      ! allocatable
      real(real_kind), allocatable :: f(:), angle(:)
      ! integer, allocatable :: work_type(:)
      ! real(real_kind), allocatable :: work_xyz(:,:), work_occ(:)

      ! NOTE: A modern optimizer should preclude the need to avoid conditionals
      ! in inner loops. It turns out that Intel vectorization routines can fail
      ! rather easily, so be sure to check the vectorization reports. Intel said
      ! they were re - working some of the vectorization to avoid this problem.
#if 0
      if (present(atom_selected)) then
         work_num_atoms = count(atom_selected/=0)
         allocate(work_xyz(3,work_num_atoms),work_type(work_num_atoms), &
               stat=alloc_status)
         REQUIRE(alloc_status==0)
         if(present(occupancy)) then
            allocate(work_occ(work_num_atoms),stat=alloc_status)
            REQUIRE(alloc_status==0)
         end if

         i = 0; iatom = 1
         do while(i<num_atoms)
            if (atom_selected(iatom)/=0) then
               i = i + 1
               work_xyz(:,i) = xyz(:,iatom)
               work_type(i) = scatter_type_index(iatom)
               if(present(occupancy)) work_occ(i) = occupancy(iatom)
            end if
            iatom = iatom + 1
         end do
         if (present(occupancy)) then
            call fourier_fcalc_no_occupancy(work_num_atoms,work_xyz,work_type)
         else
            call fourier_fcalc_with_occupancy(work_num_atoms,work_xyz,work_type,work_occ)
         end if
      else
#endif
         if (present(occupancy)) then
            call fourier_fcalc_no_occupancy(num_atoms,xyz,scatter_type_index)
         else
            call fourier_fcalc_with_occupancy(num_atoms,xyz,scatter_type_index,occupancy)
         end if
#if 0
      end if
#endif

      allocate(angle(num_atoms), f(num_atoms), stat=alloc_status)
      REQUIRE(alloc_status==0)
      return
   contains
      subroutine fourier_fcalc_with_occupancy(natom,xyz,type,occupancy)
         integer, intent(in) :: natom
         real(real_kind), intent(in) :: xyz(3,natom)
         integer, intent(in) :: type(natom)
         real(real_kind), intent(in) :: occupancy(natom)
         do ihkl = 1, num_hkl
            if (present(hkl_selected)) then
               if (hkl_selected(ihkl)==0) cycle
            end if
            ! NOTE: the atomic scatter factors do not change for a given atom type and
            ! hkl index as long as the unit cell is unchanged. If the number of atom types
            ! is small, it may be worth saving a pre - calculated array of num_type by num_reflections.
            do i = 1, num_scatter_types
               atomic_scatter_factor(i) = &
                     atom_scatter_factor_mss4(scatter_coefficients(:,:,i),mSS4(ihkl))
            end do

            ! Fhkl = SUM( fj * exp(2 * M_PI * i * (h * xj + k * yj + l * zj)) ), j = 1,num_selected_atoms
            ! where:
            !    The sum is versus j, over all selected atoms
            !    fj is the atomiuc scatter for atom j:   atomic_scatter_factor(j)
            !    h,k,l are from the hkl list for index ihkl:   hkl(1:3,ihkl)
            !    x,y,z are coordinates for atom j:   xyz(1:3,j)
            !        xyz(:) may be the actual coordinate array argument, or a reduced list.
            !
            ! Rather than using a complex exponential where the real part is always zero, this
            ! is optimized to calculate sin and cosine parts, then convert to a complex number
            ! after the A and B components are summed over all selected atoms.
            ! This can be written as:
            !
            ! Ahkl = SUM( fj * cos(2 * M_PI * (h * xj + k * yj + l * zj)) ), j = 1,num_selected_atoms
            ! Bhkl = SUM( fj * sin(2 * M_PI * (h * xj + k * yj + l * zj)) ), j = 1,num_selected_atoms

            f(:) = mSS4(ihkl) * tempFactor(:) * atomic_scatter_factor(type(:)) * occupancy(:)
            angle(:) = matmul(M_TWOPI * hkl(1:3,ihkl),xyz(1:3,:))
            Fcalc(ihkl) = cmplx( sum(f(:) * sin(angle)), sum(f(:) * cos(angle(:))) )
         end do
      end subroutine fourier_fcalc_with_occupancy

      subroutine fourier_fcalc_no_occupancy(natom,xyz,type)
         integer, intent(in) :: natom
         real(real_kind), intent(in) :: xyz(3,natom)
         integer, intent(in) :: type(natom)
         do ihkl = 1, num_hkl
            if (present(hkl_selected)) then
               if (hkl_selected(ihkl)==0) cycle
            end if
            do i = 1, num_scatter_types
               atomic_scatter_factor(i) = atom_scatter_factor_mss4(scatter_coefficients(:,:,i),mSS4(ihkl))
            end do
            f(:) = mSS4(ihkl) * tempFactor(:) * atomic_scatter_factor(type(:))
            angle(:) = matmul(M_TWOPI * hkl(1:3,ihkl),xyz(1:3,:))
            Fcalc(ihkl) = cmplx( sum(f(:) * sin(angle)), sum(f(:) * cos(angle(:))) )
         end do
      end subroutine fourier_fcalc_no_occupancy
   end subroutine fourier_Fcalc

   ! -------------------------------------------------------------------------
   ! Calculate dXYZ, dTempFactor and/or dOccupancy with respect to dF,
   ! using the direct Fourier sum.

   ! Note Q is short for occupancy; B is short for tempFactor (aka B - factor).
   ! Small molecule software uses a U parameter in place of B - factor.
   ! U has a physical meaning, but B - factor is a more natural fit to Fouriers.

   !     isotropic B - factor = 8 * pi ** 2 * isotropic - U
   !     isotropic U = [U(1,1) + U(2,2) + U(3,3)]/3.0

   subroutine fourier_dXYZBQ_dF( &
            num_hkl,hkl,dF,mSS4,hkl_selected, & ! reflections
            num_atoms,xyz,tempFactor,occupancy,scatter_type_index, & ! atoms
            num_coeffs, coeffs, & ! scatter coeff.
            dxyz,d_occupancy,d_tempFactor & ! output derivatives
      )   ! TODO: d_aniso_Bij  (anisotropic B refinement)
      implicit none
      ! reciprocal space arrays:
      integer, intent(in) :: num_hkl
      integer, intent(in) :: hkl(3,num_hkl)
      complex(real_kind), intent(in) :: dF(num_hkl)
      real(real_kind), intent(in) :: mSS4(num_hkl)
      integer, intent(in), optional :: hkl_selected(num_hkl)

      ! coordinate arrays:
      integer, intent(in) :: num_atoms
      real(real_kind), intent(in) :: xyz(3,num_atoms)
      real(real_kind), intent(in) :: tempFactor(num_atoms)
      real(real_kind), intent(in), optional :: occupancy(num_atoms)
      integer, intent(in) :: scatter_type_index(num_atoms)
      real(real_kind), intent(out), optional :: dxyz(3,num_atoms)
      real(real_kind), intent(out), optional :: d_occupancy(num_atoms)
      real(real_kind), intent(out), optional :: d_tempFactor(num_atoms)
      !real(real_kind), intent(out), optional :: d_aniso_Bij(6,num_atoms)

      ! scatter coefficient array:
      integer, intent(in) :: num_coeffs
      real(real_kind), intent(in) :: coeffs(2,scatter_ncoeffs,num_coeffs)

      ! locals
      integer :: ihkl, iatom, i
      real(real_kind), allocatable :: atomic_scatter_factor(:)
      real(real_kind) :: dhkl(3)
      !complex(real_kind) :: f(size(hkl,1))
      complex(real_kind) :: f
      real(real_kind) :: phase

      ! This array holds the scatter factor for each atom at a given HKL index.
      allocate(atomic_scatter_factor(num_scatter_types))

      REFLECTION: do ihkl = 1,num_hkl
         if (present(hkl_selected)) then
            if (hkl_selected(ihkl)==0) cycle REFLECTION
         end if

         do i = 1,num_scatter_types
            atomic_scatter_factor(i) = &
                  atom_scatter_factor_mss4(coeffs(:,:,i),mSS4(ihkl))
         end do

         ! FIXME: symmetry operations are included here, and require an additional loop.
         ! This code is currently limited to P1.
         dhkl = hkl(:,ihkl) * M_TWOPI ! * symmop...

         ! ----------------------------------------------------------------------------
         ATOM: do iatom = 1,num_atoms
            phase = sum( dhkl * xyz(:,iatom) )
            f = atomic_scatter_factor(scatter_type_index(iatom)) &
                  * exp(mSS4(ihkl) * tempFactor(iatom)) &
                  * cmplx(cos(phase),sin(phase))

            if (present(d_occupancy)) then
               d_occupancy(iatom) = d_occupancy(iatom) + real(f) * real(dF(ihkl)) + aimag(f) * aimag(dF(ihkl))
            end if

            if (present(occupancy)) then
               f = f * occupancy(iatom)
            end if

            if (present(d_tempFactor)) then
               d_tempFactor(iatom) = d_tempFactor(iatom) &
                     + ( real(f) * real(dF(ihkl)) + aimag(f) * aimag(dF(ihkl)) ) &
                     * mSS4(ihkl)
            end if

            if (present(dxyz)) then
               dxyz(:,iatom) = dxyz(:,iatom) + dhkl * &
                     ( real(f) * aimag(dF(ihkl)) - aimag(f) * real(dF(ihkl)) )
            end if

         end do ATOM

      end do REFLECTION

   end subroutine fourier_dXYZBQ_dF

   ! -----------------------------------------------------------------------------
   ! R - factor = sum(abs(Fobs - Fcalc)) / sum(Fobs)
   ! This routine computes the force gradient on Fcalc as a harmonic
   ! restraint
   subroutine dF_dTarget(fobs,Fcalc,weight,selected,deriv,residual)
      implicit none
      real(real_kind), intent(in) :: fobs(:)
      complex(real_kind), intent(in) :: Fcalc(:)
      integer, intent(in), optional :: selected(:)
      real(real_kind), intent(in), optional :: weight (:)
      complex(real_kind), intent(out), optional :: deriv(:)
      real(real_kind), intent(out), optional :: residual

      integer :: num_hkl
      real(real_kind) :: sum_fo_fc, sum_fo_squared, sum_fc_squared
      real(real_kind), allocatable :: abs_Fcalc(:)
      real(real_kind) :: Fcalc_scale, norm_scale
      real(real_kind), parameter :: F_EPSILON = 1.0e-20_rk_

      num_hkl = size(fobs,1)
      allocate(abs_Fcalc(num_hkl))
      abs_Fcalc = abs(Fcalc)

      if (present(weight)) then
         if (present(selected)) then
            sum_fo_fc = sum(weight * fobs * abs_Fcalc,selected/=0)
            sum_fo_squared = sum(weight * fobs ** 2,selected/=0)
            sum_fc_squared = sum(weight * abs_Fcalc ** 2,selected/=0)
         else
            sum_fo_fc = sum(weight * fobs * abs_Fcalc)
            sum_fo_squared = sum(weight * fobs ** 2)
            sum_fc_squared = sum(weight * abs_Fcalc ** 2)
         end if
      else
         if (present(selected)) then
            sum_fo_fc = sum(fobs * abs_Fcalc,selected/=0)
            sum_fo_squared = sum(fobs ** 2,selected/=0)
            sum_fc_squared = sum(abs_Fcalc ** 2,selected/=0)
         else
            sum_fo_fc = sum(fobs * abs_Fcalc)
            sum_fo_squared = sum(fobs ** 2)
            sum_fc_squared = sum(abs_Fcalc ** 2)
         end if
      end if

      if (sum_fc_squared<F_EPSILON .or. sum_fo_squared<F_EPSILON) then
         if (present(deriv)) then
            if (present(selected)) then
               where (selected/=0) deriv = (0.0_rk_,0.0_rk_)
            else
               deriv = (0.0_rk_,0.0_rk_)
            end if
         end if
         if (present(residual)) residual = -1.0_rk_
         return
      end if

      Fcalc_scale = sum_fo_fc / sum_fc_squared
      norm_scale = 1.0_rk_ / sum_fo_squared

      !write( * ,*) 'SUMS = ',sum_fo_fc, sum_fo_squared, sum_fc_squared
      !write( * ,*) 'scale = ',norm_scale,', Fcalc_scale = ',Fcalc_scale

      ! Note: when Fcalc is approximately zero the phase is undefined, so no force
      ! can be determined even if the energy is high. (Similar to a linear bond angle.)

      ! We get the phase from Fcalc, so it needs to be divided by abs(Fcalc)
      ! to become a unit vector.
      !
      ! deriv = Fcalc/abs(Fcalc) * K * ( Fobs - Fcalc_scale * abs(Fcalc) )
      !
      !       = Fcalc * K * ( Fobs/abs(Fcalc) - Fcalc_scale )

      if (present(deriv)) then
         if (present(weight)) then
            if (present(selected)) then
               where (selected/=0)
               where (abs_Fcalc > F_EPSILON)
               deriv = Fcalc * ( - 2.0_rk_ * Fcalc_scale * &
                     weight * norm_scale * ( fobs / abs_Fcalc - Fcalc_scale ) )
            else where
               deriv = (0.0_rk_,0.0_rk_)
               end where
               end where
            else ! no selected
               where (abs_Fcalc > F_EPSILON)
               deriv = Fcalc * ( - 2.0_rk_ * Fcalc_scale * &
                     weight * norm_scale * ( fobs / abs_Fcalc - Fcalc_scale ) )
            else where
               deriv = (0.0_rk_,0.0_rk_)
               end where
            end if
         else ! no weight
            if (present(selected)) then
               where (selected/=0)
               where (abs_Fcalc > F_EPSILON)
               deriv = Fcalc * ( - 2.0_rk_ * Fcalc_scale * &
                     norm_scale * ( fobs/abs_Fcalc - Fcalc_scale ) )
            else where
               deriv = (0.0_rk_,0.0_rk_)
               end where
               end where
            else ! no selected (and no weight)
               where (abs_Fcalc > F_EPSILON)
               deriv = Fcalc/abs_Fcalc * ( - 2.0_rk_ * Fcalc_scale * &
                     norm_scale * ( fobs - Fcalc_scale * abs_Fcalc ) )
               !deriv = Fcalc * ( - 2.0_rk_ * Fcalc_scale * &
                     !      norm_scale * ( fobs / abs_Fcalc - Fcalc_scale ) )
            else where
               deriv = (0.0_rk_,0.0_rk_)
               end where
            end if
         end if
      end if

      ! Do the appropriate R - factor calculation, depending on weights and/or a selection mask.
      if (present(residual)) then
         if (present(weight)) then
            if (present(selected)) then
               residual = sum(weight * abs(fobs - Fcalc_scale * abs_Fcalc),selected/=0) &
                     / sum(fobs,selected/=0)
            else
               residual = sum (weight * abs( fobs - Fcalc_scale * abs_Fcalc )) / sum(fobs)
            end if
         else
            if (present(selected)) then
               residual = sum(abs(fobs - Fcalc_scale * abs_Fcalc), selected/=0) / sum(fobs, selected/=0)
               ! TARGET ENERGY = residual = norm_scale &
                     !          * sum((fobs - Fcalc_scale * abs_Fcalc) ** 2, selected/=0) / count(selected/=0)
            else
               residual = sum (abs(fobs - Fcalc_scale * abs_Fcalc)) / sum(fobs)
            end if
         end if
      end if

   end subroutine dF_dTarget

   ! ==========================================================================

   subroutine FFT_Fcalc( &
            !num_hkl,hkl,Fcalc,mSS4,hkl_selected, & ! reflections
            num_hkl,Fcalc,hkl_selected, & ! reflections
      num_atoms,xyz,tempFactor,occupancy,scatter_type, & ! atoms
      num_types,num_coeff,coeffs) ! scatter coeff.
      use xray_utils_module, only: allocate_lun, write_map_ezd
      implicit none

      integer, intent(in) :: num_hkl
      ! integer, intent(in) :: hkl(3,num_hkl)                ! H,K,L indices
      complex(real_kind), intent(out) :: Fcalc(num_hkl)    ! Fcalc(h,k,l)
      ! real(real_kind), intent(in) :: mSS4(num_hkl)         ! -(S^2/4)
      integer, intent(in), optional :: hkl_selected(num_hkl)

      integer, intent(in) :: num_atoms
      real(real_kind), intent(in) :: xyz(3,num_atoms) !fractional
      real(real_kind), intent(in) :: tempFactor(num_atoms)
      real(real_kind), intent(in), optional :: occupancy(num_atoms)
      integer, intent(in) :: scatter_type(num_atoms)

      integer, intent(in) :: num_coeff, num_types
      real(real_kind), intent(in) :: coeffs(2,num_coeff,num_types)

      ! integer :: i, j, ih,ik,il, unit
      integer :: unit

      ! Calculate the high-res limit
      !FIXME: this calculation is wrong.
      !if (present(hkl_selected)) then
      !  high_res_limit = 1.0_rk_/sqrt( -4.0 * minval(mSS4,mask = hkl_selected))
      !else
      !  high_res_limit = 1.0_rk_/sqrt( -4.0 * minval(mSS4))
      !end if

      density_map = 0

      call grid_atomic_density( &
            num_atoms, xyz, tempFactor, occupancy, scatter_type, &
            num_types, num_coeff, coeffs)

      call amopen(allocate_lun(unit),'xray_atom.map','U','F','R')
      call write_map_ezd(unit, &
          reshape((/0,grid_size(1)-1,0,grid_size(2)-1,0,grid_size(3)/),(/2,3/)), &
          density_map)
      stop

      !FIXME: call FFT_forward(FFT)

      call FFT_grid_to_hkl(density_map,num_hkl,hkl_index,Fcalc,hkl_selected)
#if 0
      ! Extract reflections from the FFT that are present in our selected HKL indices.
      do i = 1,num_hkl
         if (present(hkl_selected)) then
            if (hkl_selected(i)==0) cycle
         end if

         ! The FFT calculates the zero - plane and the half where H > 0
         ! Use point symmetry to ensure H >= 0
         ! Point inversion requires converting to the complex conjugate.
         ! NOTE: From testing, the need for complex-conjugate is reversed from
         ! what I expected. (TODO: Look over the math again.)
         ! Negative K,L values wrap at the origin.
         if (hkl(1,i)<0) then
            ih = -hkl(1,i)
            ik = -hkl(2,i); if (ik<0) ik = ik + recip_size(2)
            il = -hkl(3,i); if (il<0) il = il + recip_size(3)
            Fcalc(i) = fft_normalization * cmplx(recip_grid(1,ih,ik,il),recip_grid(1,ih,ik,il)) &
                  * exp( fft_bfactor_sharpen * mSS4(i))
         else
            ih = hkl(1,i)
            ik = hkl(2,i); if (ik<0) ik = ik + recip_size(2)
            il = hkl(3,i); if (il<0) il = il + recip_size(3)
            ! complex-conjugate:
            Fcalc(i) = fft_normalization * cmplx(recip_grid(1,ih,ik,il),-recip_grid(1,ih,ik,il)) &
                  * exp( fft_bfactor_sharpen * mSS4(i))
         end if
      end do
#endif
   end subroutine FFT_Fcalc

   subroutine FFT_dXYZBQ_dF(&
            !num_hkl,hkl,dF,mSS4,hkl_selected, & ! reflections
            num_hkl,dF,hkl_selected, & ! reflections
      num_atoms,xyz,tempFactor,occupancy,scatter_type, & ! atoms
      num_coeff, num_types, coeffs, & ! scatter coeff.
      dxyz,d_occupancy,d_tempFactor & ! output derivatives
      )   ! TODO: d_aniso_Bij

      implicit none

      ! reciprocal space arrays:
      integer, intent(in) :: num_hkl
      ! integer, intent(in) :: hkl(3,num_hkl)
      complex(real_kind), intent(in) :: dF(num_hkl)
      ! real(real_kind), intent(in) :: mSS4(num_hkl)
      integer, intent(in), optional :: hkl_selected(num_hkl)

      ! coordinate arrays:
      integer, intent(in) :: num_atoms
      real(real_kind), intent(in) :: xyz(3,num_atoms)
      real(real_kind), intent(in) :: tempFactor(num_atoms)
      real(real_kind), intent(in), optional :: occupancy(num_atoms)
      integer, intent(in) :: scatter_type(num_atoms)
      real(real_kind), intent(out), optional :: dxyz(3,num_atoms)
      real(real_kind), intent(out), optional :: d_occupancy(num_atoms)
      real(real_kind), intent(out), optional :: d_tempFactor(num_atoms)
      !real(real_kind), intent(out), optional :: d_aniso_Bij(6,num_atoms)

      ! scatter coefficient array:
      integer, intent(in) :: num_coeff, num_types
      real(real_kind), intent(in) :: coeffs(num_coeff,2,num_types)

      ! integer :: ihkl, isym, i, j, ih,ik,il
      ! real(real_kind) :: scale

      ! This array holds the scatter factor for each atom at a given HKL index.
      ! Alternatively, atoms could be sorted by type, but that is likely inefficient,
      ! because sorting is slow and that order would be inefficient for non - bond forces.
      ! -----------------------------------------------------------------------------
      ! Map reflections to the grid, applying symmetry operations:
      call FFT_hkl_to_grid(density_map,num_hkl,hkl_index,dF,hkl_selected)
#if 0
      scale = 1.0_rk_/fft_normalization
      recip_grid = (0.0_rk_,0.0_rk_)
      do isym = 1,num_symmops
         do ihkl = 1,num_hkl
            if (present(hkl_selected)) then
               if (hkl_selected(ihkl)==0) cycle
            end if
            ! Note: H==0 plane requires both +/- symmetry points.
            if (hkl(1,ihkl) <= 0) then
               ih = -hkl(1,ihkl)
               ik = -hkl(2,ihkl); if (ik<0) ik = ik + recip_size(2)
               il = -hkl(3,ihkl); if (il<0) il = il + recip_size(3)
               recip_grid(ih,ik,il) = recip_grid(ih,ik,il) &
                     + dF(ihkl) * scale / exp( fft_bfactor_sharpen * mSS4(ihkl))
            end if
            if (hkl(1,ihkl) >= 0) then
               ik = hkl(1,i)
               ik = hkl(2,i); if (ik<0) ik = ik + recip_size(2)
               il = hkl(3,i); if (il<0) il = il + recip_size(3)
               recip_grid(ih,ik,il) = recip_grid(ih,ik,il) + conjg(dF(ihkl) &
                     * scale / exp( fft_bfactor_sharpen * mSS4(ihkl)))
            end if
         end do
      end do
#endif
      !FIXME: call FFT_backward(FFT)

      call grid_atomic_density( &
            num_atoms, xyz, tempFactor, occupancy, scatter_type, &
            num_types, num_coeff, coeffs, &
            d_tempFactor = d_tempFactor, &
            d_occupancy = d_occupancy, &
            d_XYZ = dxyz)

   end subroutine FFT_dXYZBQ_dF

   ! -------------------------------------------------
   ! Converts from HKL reflection list to FFT grid.
   ! Applies symmetry operators.
   subroutine FFT_hkl_to_grid(recip_grid, &
            num_hkl,hkl_index,refl,hkl_selected)

      real(fft_kind) :: &
            recip_grid(2,0:recip_size(1)-1,0:recip_size(2)-1,0:recip_size(3-1))
      integer, intent(in) :: num_hkl
      integer, intent(in) :: hkl_index(3,num_hkl)
      complex(real_kind), intent(in) :: refl(num_hkl)
      integer, intent(in), optional :: hkl_selected(num_hkl)

      real(real_kind) :: scale
      real(real_kind) :: grid_to_phase(3), phase
      integer :: isym, ihkl, h,k,l,hm,km,lm
      real(real_kind) :: rsymm_hkl(3), s(3,3), t(3)
      integer :: symm_hkl(3)
      real(real_kind) :: r
      logical :: Friedel
      complex(real_kind) :: f

      grid_to_phase = M_TWOPI/real(grid_size,real_kind)
      scale = 1.0_rk_/fft_normalization

      recip_grid = 0.0_rk_

      ! loop over all symmetry operators
      do isym = 1,num_symmops
         ! fractional translation converted to phase-shift:
         t = symmop(1:3,4,isym) * M_TWOPI
         ! Reciprocal symmop is the transpose of the real-space symmop-op
         s = transpose(symmop(1:3,1:3,isym))

         do ihkl = 1,num_hkl
            if(present(hkl_selected)) then
               if (hkl_selected(ihkl)==0) cycle
            end if

            rsymm_hkl = matmul(real(hkl_index(1:3,ihkl)),s)
            phase = sum(rsymm_hkl * t) ! FIXME
            f = cmplx(cos(phase),sin(phase))
            ! Convert to array indices (when needed):
            symm_hkl = modulo(nint(rsymm_hkl),grid_size)

            ! Map h,k,l to the H>=0 hemisphere:
            Friedel = (symm_hkl(1) > grid_size(1)/2)
            if (Friedel) symm_hkl = modulo(-symm_hkl,grid_size)

            ! Note: the derivative of a constant factor is conjugate,
            ! as are the factors for the subsequent fourier transform.
            f = f * conjg(refl(ihkl)) &
                  * scale / exp( fft_bfactor_sharpen * mSS4(ihkl))
            if (Friedel) f = conjg(f)
            recip_grid(:,symm_hkl(1),symm_hkl(2),symm_hkl(3)) = &
                  recip_grid(:,symm_hkl(1),symm_hkl(2),symm_hkl(3)) &
                  + (/real(f),aimag(f)/)

         end do !ihkl
      end do !isym

      ! The H == 0 plane contains part of the symmetric hemisphere.
      ! The data must be symmetrically summed and distributed.
      ! From CNS: Note the factor 1/2 in the calling routine.
      do l = 0,grid_size(3)/2
         lm = modulo(l,grid_size(3))
         if (l /= lm) cycle
         do k = 0,grid_size(2)/2
            km = modulo(k,grid_size(2))

            ! add complementary conjugates of 0,k,l and 0,km,lm
            r = recip_grid(1,0,k,l) + recip_grid(1,0,km,lm)
            recip_grid(1,0,k,l) = r
            recip_grid(1,0,km,lm) = r

            r = recip_grid(2,0,k,l) - recip_grid(2,0,km,lm)
            recip_grid(1,0,k,l) = r
            recip_grid(1,0,km,lm) = -r

            do h = 1,merge(grid_size(1)/2, grid_size(1)-1, k == km)
               !recip_grid(h,k,l) = recip_grid(h,k,l) &
               !                  + conjg(recip_grid(grid_size(1)-h,km,lm))
               !recip_grid(grid_size(1)-h,km,lm) = dconjg(recip_grid(h,k,l))
               hm = grid_size(1)-h

               r = recip_grid(1,h,k,l) + recip_grid(1,hm,km,lm)
               recip_grid(1,h,k,l) = r
               recip_grid(1,hm,km,lm) = r

               r = recip_grid(2,h,k,l) - recip_grid(2,hm,km,lm)
               recip_grid(1,h,k,l) = r
               recip_grid(1,hm,km,lm) = -r

            end do
         end do
      end do
   end subroutine FFT_hkl_to_grid

   ! -------------------------------------------------
   ! Converts between FFT grid format and reflection list format.
   ! Applies symmetry operators.
   subroutine FFT_grid_to_hkl(recip_grid, &
            num_hkl,hkl_index,refl,hkl_selected)
      implicit none
      real(fft_kind) :: &
            recip_grid(2,0:recip_size(1)-1,0:recip_size(2)-1,0:recip_size(3-1))
      integer, intent(in) :: num_hkl
      integer, intent(in) :: hkl_index(3,num_hkl)
      complex(real_kind), intent(out) :: refl(num_hkl)
      integer, intent(in), optional :: hkl_selected(num_hkl)

      real(real_kind) :: grid_to_phase(3), phase
      integer :: isym, ihkl
      real(real_kind) :: rsymm_hkl(3), s(3,3), t(3)
      integer :: symm_hkl(3)
      logical :: Friedel
      complex(real_kind) :: f
      real(real_kind) :: r(2)

      grid_to_phase = M_TWOPI/real(grid_size,real_kind)

      ! fill the reflection array with initial indentity-symmop values:
      do ihkl = 1,num_hkl
         if(present(hkl_selected)) then
            if (hkl_selected(ihkl)==0) cycle ! set unselected to zero?
         end if

         rsymm_hkl = matmul(real(hkl_index(1:3,ihkl)),s)
         symm_hkl = modulo(nint(rsymm_hkl),grid_size)

         ! Map h,k,l to the H >= 0 hemisphere:
         Friedel = (symm_hkl(1) > grid_size(1)/2)
         if (Friedel) symm_hkl = modulo(-symm_hkl,grid_size)

         r = recip_grid(:,symm_hkl(1),symm_hkl(2),symm_hkl(3))
         if (Friedel) r(2) = -r(2) !conjg(r)
         refl(ihkl) = fft_normalization * cmplx(r(1),r(2))

      end do !ihkl

      ! loop over all other symmetry operators
      do isym = 2,num_symmops
         t = symmop(1:3,4,isym) * M_TWOPI
         s = transpose(symmop(1:3,1:3,isym))

         do ihkl = 1,num_hkl
            if(present(hkl_selected)) then
               if (hkl_selected(ihkl)==0) cycle
            end if

            ! compute the phase shift from the reciprocal
            ! symmetry index (always zero for P1)
            rsymm_hkl = matmul(real(hkl_index(1:3,ihkl)),s)
            phase = sum(rsymm_hkl * t) ! FIXME!
            f = cmplx(cos(phase),sin(phase))

            ! Convert to array indices (when needed):
            symm_hkl = modulo(nint(rsymm_hkl),grid_size)

            ! Map h,k,l to the H >= 0 hemisphere:
            Friedel = (symm_hkl(1) > grid_size(1)/2)
            if (Friedel) symm_hkl = modulo(-symm_hkl,grid_size)

            r = recip_grid(:,symm_hkl(1),symm_hkl(2),symm_hkl(3))
            if (Friedel) r(2) = -r(2) !conjg(r)
            refl(ihkl) = refl(ihkl) + f * fft_normalization * cmplx(r(1),r(2))

         end do !ihkl
      end do !isym

   end subroutine FFT_grid_to_hkl

   function atom_scatter_factor_mss4(coeffs,mss4) result(sfac)
      real(real_kind) :: sfac
      real(real_kind), intent(in) :: coeffs(2,scatter_ncoeffs), mss4
      sfac = coeffs(1,scatter_ncoeffs) + &
             sum( coeffs(1,1:scatter_ncoeffs-1) &
             * exp(mss4*coeffs(2,1:scatter_ncoeffs-1)))
   end function atom_scatter_factor_mss4

   !=====================================================================
   ! Increment NUM to the next value with no factors greater than
   ! FFT_MAX_PRIME, that also contains the factor FACTOR.
   pure subroutine fft_adjust_size (num,factor)
      implicit none
      integer, intent(inout) :: num
      integer, intent(in), optional :: factor ! required factor
      do
         if (mod(num,factor)==0) then
            if (max_factor(num,fft_max_prime)) return
         end if
         num = num + 1
      end do
   end subroutine fft_adjust_size

   !=====================================================================
   ! Return FALSE if any factors of NUM are greater than FFT_MAX_PRIME
   pure function max_factor(num, max) result(ok)
      implicit none
      logical :: ok
      integer, intent(in) :: num, max
      ! local
      integer :: n, p
      ! begin
      ok = .true.
      if (num <= 1) return ! OK (assumes no negative numbers)
      if (max < 2) then
         ok = .false.
         return
      end if
      n = num
      do while (.not.btest(n,0))
         n = n / 2
      end do
      if (n == 1) return
      do p = 3, max, 2
         do while (mod(n, p) == 0)
            n = n / p
         end do
         if (n == 1) return
      end do
      ok = .false.
      return
   end function max_factor
   !=======================================================================

   subroutine test_fft()
      implicit none
      ! Test FFT library to ensure it actually works.
      ! local
      integer :: i, u, v, w, n, alloc_status
      logical :: error
      integer :: grid_size(3) !, rm(1)
      type(fft1c_descriptor) :: cfft
      type(fft1cr_descriptor) :: crfft
      type(fft3c_descriptor) :: cfft3
      logical, parameter :: verbose=.true.
      real(fft_kind) :: err
      ! allocatable
      complex(fft_kind), allocatable, target :: cmap(:,:,:), cmap_orig(:,:,:)
      ! complex/real equivalence
      integer, parameter :: MAXSEQ=256
      complex(fft_kind) :: cseq(0:MAXSEQ-1)
      real(fft_kind) :: rseq(0:MAXSEQ*2-1)
      equivalence(cseq,rseq)
      real(fft_kind) :: rseq_orig(0:MAXSEQ-1)

      ! parameter
      real(fft_kind), parameter :: max_err = 0.01
      real(fft_kind), parameter :: rzero(1) = (/0.0/)
      complex(fft_kind), parameter :: czero(1) = rzero

      ! begin

      ! Set up complex-to-complex transform
      grid_size = (/ 4*3, 3*5, 5*2 /)
      do i = 1, 3
         call fft_adjust_size(grid_size(i))
      end do
      allocate(cmap(grid_size(1),grid_size(2),grid_size(3)), &
            cmap_orig(grid_size(1),grid_size(2),grid_size(3)), &
            stat=alloc_status)
      REQUIRE(alloc_status==0)

      ! FFT3C -> conjugate complex -> FFT3C -> VFYCMAP
      ! Fill CMAP with a pattern:
      do u = 1, grid_size(1)
         do v = 1, grid_size(2)
            do w = 1, grid_size(3)
               cmap(u,v,w) = cmplx( &
                     real(mod(u+v+w,11),fft_kind), &
                     real(mod(u+v+w,13),fft_kind))
            end do
         end do
      end do
      cmap_orig = cmap

      call fft3c_new(cfft3,grid_size,error)
      if (error) then
         write(stderr,'(A)') 'Fatal Error: Call FFT3C_NEW failed.'
         call mexit(stderr,1)
      end if
      call fft3c(cfft3,cmap, error)
      if (error) then
         write(stderr,'(A)') 'Fatal Error: Call FFT3C 1 failed.'
         call mexit(stderr,1)
      end if
      cmap=conjg(cmap)/real(product(grid_size),fft_kind)
      call fft3c(cfft3,cmap, error)
      if (error) then
         write(stderr,'(A)') 'Fatal Error: Call FFT3C 2 failed.'
         call mexit(stderr,1)
      end if

      ! See if we got the right numbers back
      err = max(maxval(abs(real(cmap)-real(cmap_orig))), &
            maxval(abs(aimag(cmap)-aimag(cmap_orig))))
      if (err > max_err) then
         write(stderr,'(A)') 'Fatal Error: Test of FFT3C failed.'
         call mexit(stderr,1)
      end if
      deallocate(cmap,cmap_orig)

      ! Test single precision complex-to-complex-to-real transform
      n = 1
      do
         n = n + 1
         call fft_adjust_size(n,2)
         if (n>MAXSEQ) exit

         !call mfftcr(1, nn, rm)
         ! SFFT1C -> SFFT1CR (with unique half) -> VFYSSEQ
         !!!call isseq(n,cseq)
         do i = 0,n-1
            rseq_orig(i) = real(mod(i+4,7),fft_kind)
         end do
         cseq = rseq_orig

         call fft1c_new(cfft,n,error)
         if (error) then
            write(stderr,'(A)') 'Fatal Error: CALL SFFT1C 1 failed.'
            call mexit(6,1)
         end if
         call fft1c(cfft, cseq)
         call fft1c_delete(cfft)

         call fft1cr_new(crfft,n,error)
         if (error) then
            write(stderr,'(A)') 'Fatal Error: CALL SFFT1CR 1 failed.'
            call mexit(6,1)
         end if
         call fft1cr(crfft, rseq)
         call fft1cr_delete(crfft)

         ! transformed real is scaled and reversed.
         error = any(nint(rseq(n-1:0)/real(n,fft_kind))/=nint(rseq_orig))
         if (error) then
            write(stderr, '(1X,A,I3)') &
                  'Fatal Error: Test of SFFT1C and SFFT1CR failed for N = ',n
            exit
         end if
         if (verbose) then
            write(stdout, '(1X,A,I3)') &
                  'Test of SFFT1C and SFFT1CR succeeded for N = ',n
         end if
      end do
      return
   end subroutine test_fft

end module xray_fourier_module