multidip_smooth_rawdips.f90

Go to the documentation of this file.

! Copyright 2020
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-out (UKRmol+ suite).
!
!     UKRmol-out is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
!
!     UKRmol-out is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with  UKRmol-out (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
!
 
!> \brief   Smoothing of MULTIDIP raw transition dipoles
!> \author J Benda
!> \date    2021
!>
!> Program usage:
!>
!>     smooth_rawdips rawdip-*.txt
!>
!> Applies Gaussian smoothing to the raw dipole elements written by MULTIDIP, producing "XXX-smooth.txt" file for every "XXX.txt"
!> given on command line. The smoothing bandwidth is proportional to the global parameter `sigma`. The smoothing is performed in
!> iterations, starting with identically zero dataset. In each iteration, every sample is assigned a weight that is approximately
!> inversely proportional to its distance from the current smoothed value at the same energy. This enables the program to produce
!> results that are not virtually unaffected by extremely thin but extremely large spikes, which would otherwise affect the
!> Gaussian mean in a significant way. The number of these smoothing iterations is controlled by the global parameter `niter`.
!>
program multidip_smooth_rawdips
 
    use iso_fortran_env, only: real64
 
    implicit none
 
    real(real64), parameter :: sigma = 1e-3     ! Gaussian filter sharpness (large values will conserve narrow features)
    integer,      parameter :: niter = 5        ! number of smoothing iterations (to compensate large out-of-trend spikes)
 
    real(real64), allocatable :: dips(:, :)
    character(len=1024)       :: rawname, smoothname
    integer                   :: iarg, narg, dot
 
    narg = command_argument_count()
 
    if (narg == 0) then
        print '(A)', 'Just give me a list of files to smooth...'
        stop
    end if
 
    do iarg = 1, narg
 
        ! get next filename from the command line
        call get_command_argument(iarg, rawname)
        smoothname = rawname
 
        ! translate 'XXX.YYY' to 'XXX-smooth.YYY'
        dot = scan(rawname, '.', back = .true.)
        if (dot /= 0) then
            smoothname = smoothname(1:dot-1) // '-smooth' // smoothname(dot:)
        end if
 
        print '(3A)', trim(rawname), ' -> ', trim(smoothname)
 
        ! perform the smoothing
        call read_dipoles(trim(rawname), dips)
        call smooth_dipoles(dips)
        call write_dipoles(trim(smoothname), dips)
 
    end do
 
contains
 
    subroutine read_dipoles (filename, dips)
 
        character(len=*),          intent(in)    :: filename
        real(real64), allocatable, intent(inout) :: dips(:, :)
 
        integer :: i, u, ierr, n
        real(real64) :: x, y
 
        open (newunit = u, file = filename, action = 'read', form = 'formatted')
 
        ! first, count the lines
        n = 0
        ierr = 0
        do while (ierr == 0)
            read (u, *, iostat = ierr) x, y
            if (ierr == 0) then
                n = n + 1
            end if
        end do
 
        ! set up dipole storage
        if (allocated(dips)) then
            if (size(dips, 2) /= n) then
                deallocate (dips)
            end if
        end if
        if (.not. allocated(dips)) then
            allocate (dips(2, n))
        end if
 
        ! now read from beginning
        rewind(u)
        do i = 1, n
            read (u, *) dips(:, i)
        end do
 
        close (u)
 
        ! replace NaNs with huge values, which will have negligible effect on the smoothing
        where (.not. dips == dips)
            dips = huge(dips)
        end where
 
    end subroutine read_dipoles
 
 
    subroutine write_dipoles (filename, dips)
 
        character(len=*),          intent(in) :: filename
        real(real64), allocatable, intent(in) :: dips(:, :)
 
        integer :: u
 
        open (newunit = u, file = filename, action = 'write', form = 'formatted')
 
        write (u, '(2E25.15E3)') dips
 
        close (u)
 
    end subroutine write_dipoles
 
 
    subroutine smooth_dipoles (dips)
 
        real(real64), allocatable, intent(inout) :: dips(:, :)
        real(real64), allocatable :: weights(:), dips0(:, :), dips1(:, :)
        real(real64) :: weight, weight_sum, re_d, im_d
        integer :: i, j, iter, n0, n, d
 
        n = size(dips, 2)
 
        allocate (weights(n), dips0(2, n), dips1(2, n))
 
        ! set up the Gaussian weights
        do i = 1, n
            weights(i) = exp(-sigma*(i - 1)**2)
            dips0(:, i) = 0
        end do
 
        ! count the number of leading zeros; they will not participate in smoothing
        n0 = 1
        do i = 1, n
            if (any(dips(:, i) /= 0)) then
                n0 = i
                exit
            end if
        end do
 
        ! do the smoothing iterations
        do iter = 1, niter
            ! calculate entries of the smoothed dataset
            !$omp parallel do default(none) private(i,j,weight_sum,d,re_d,im_d,weight) shared(dips,dips1,dips0,weights,n,n0)
            do i = n0, n
                weight_sum = 0
                dips1(:, i) = 0
                ! use all elements from the previous iteration
                do j = n0, n
                    d = 1 + abs(i - j)
                    re_d = dips(1, j) - dips0(1, j)
                    im_d = dips(2, j) - dips0(2, j)
                    weight = weights(d) / sqrt(1 + re_d*re_d + im_d*im_d)
                    dips1(:, i) = dips1(:, i) + dips(:, j) * weight
                    weight_sum = weight_sum + weight
                end do
                dips1(:, i) = dips1(:, i) / weight_sum
            end do
            !$omp end parallel do
            dips0 = dips1
        end do
        dips = dips0
 
    end subroutine smooth_dipoles
 
end program multidip_smooth_rawdips