sk.F90

Go to the documentation of this file.

! This file is part of the 'atomes' software.
!
! 'atomes' is free software: you can redistribute it and/or modify it under the terms
! of the GNU Affero General Public License as published by the Free Software Foundation,
! either version 3 of the License, or (at your option) any later version.
!
! 'atomes' 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 Affero General Public License along with 'atomes'.
! If not, see <https://www.gnu.org/licenses/>
!
! Copyright (C) 2022-2026 by CNRS and University of Strasbourg
!
!>
!! @file sk.F90
!! @short S(k) analysis: direct reciprocal space calculation
!! @author Sébastien Le Roux <sebastien.leroux@ipcms.unistra.fr>
 
INTEGER (KIND=c_int) FUNCTION s_of_k (NQ, XA) bind (C,NAME='s_of_k_')
 
! Total structure factor
! Partial structure factors from k-points
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: nq, xa
DOUBLE PRECISION :: factor, xfactor
 
INTERFACE
  DOUBLE PRECISION FUNCTION fqx(TA, Q)
    INTEGER, INTENT(IN) :: ta
    DOUBLE PRECISION, INTENT(IN) :: q
  END FUNCTION
END INTERFACE
 
if(allocated(sij)) deallocate(sij)
allocate(sij(nq,nsp,nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k"//char(0), "Table: Sij"//char(0))
  s_of_k = 0
  goto 001
endif
sij(:,:,:)=0.0d0
 
if(allocated(cij)) deallocate(cij)
allocate(cij(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k"//char(0), "Table: cij"//char(0))
  s_of_k = 0
  goto 001
endif
if(allocated(sik)) deallocate(sik)
allocate(sik(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k"//char(0), "Table: sik"//char(0))
  s_of_k = 0
  goto 001
endif
 
#ifdef OPENMP
  !t0 = OMP_GET_WTIME ()
  call fourier_trans_qvect () ! Default Q-vector parallelization
  !t1 = OMP_GET_WTIME ()
  !write (*,*) "temps d’excecution QVT 2:", t1-t0
#else
  call fourier_trans_steps ()
#endif
 
if (allocated(cij)) deallocate(cij)
if (allocated(sik)) deallocate(sik)
if (allocated(qvectx)) deallocate(qvectx)
if (allocated(qvecty))deallocate(qvecty)
if (allocated(qvectz)) deallocate(qvectz)
if (allocated(modq)) deallocate(modq)
 
if(allocated(s)) deallocate(s)
allocate(s(nq), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k"//char(0), "Table: S"//char(0))
  s_of_k = 0
  goto 001
endif
if(allocated(xs)) deallocate(xs)
allocate(xs(nq), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k"//char(0), "Table: XS"//char(0))
  s_of_k = 0
  goto 001
endif
 
factor=0.0d0
do i=1, nsp
  factor=factor + nbspbs(i)*nscattl(i)**2
enddo
factor=factor*ns
 
if (xa .eq. 1) then
  xfactor=0.0d0
  do i=1, nsp
    xfactor=xfactor + nbspbs(i)*xscattl(i)**2
  enddo
  xfactor=xfactor*ns
endif
 
s(:)=0.0d0
xs(:)=0.0d0
 
do i=1,  nq
  if (degeneracy(i) .gt. 0) then
    do k=1, nsp
      do j=1, nsp
        s(i)=s(i)+sij(i,k,j)*nscattl(k)*nscattl(j)
        if (xa .eq. 1) then
          xs(i)=xs(i)+sij(i,k,j)*xscattl(k)*xscattl(j)
        else
          xs(i)=xs(i)+sij(i,k,j)*fqx(int(xscattl(k)),k_point(i))*fqx(int(xscattl(j)),k_point(i))
        endif
      enddo
    enddo
    s(i)=s(i)/(factor*degeneracy(i))
    if (xa .eq. 1) then
      xs(i)=xs(i)/(xfactor*degeneracy(i))
    else
      xfactor=0.0d0
      do k=1, nsp
        xfactor=xfactor + nbspbs(k)*fqx(int(xscattl(k)),k_point(i))**2
      enddo
      xs(i)=xs(i)/(xfactor*degeneracy(i)*ns)
    endif
    do j=1, nsp
      do k=1, nsp
        sij(i,j,k)=sij(i,j,k)/(degeneracy(i)*sqrt(dble(nbspbs(k)*nbspbs(j)))*ns)
      enddo
    enddo
  endif
enddo
 
s_of_k = sk_save()
 
001 continue
 
if (allocated(degeneracy)) deallocate(degeneracy)
if (allocated(k_point)) deallocate(k_point)
if (allocated(sij)) deallocate(sij)
if (allocated(s)) deallocate(s)
if (allocated(xs)) deallocate(xs)
 
CONTAINS
 
!************************************************************
!
! this subroutine computes the sine and cosine sums from the
! configuration for all q-vectors needed
!
SUBROUTINE fourier_trans_steps ()
 
  USE parameters
 
  IMPLICIT NONE
  DOUBLE PRECISION :: qx, qy, qz, qtr, sini, cosi
  INTEGER :: sp
 
  do k=1, ns
    do j=1, number_of_qvect
 
      l=anint((modq(j)-qvmin)/delta_q)+1
 
      if (l .le. nq) then
 
        qx=qvectx(j)
        qy=qvecty(j)
        qz=qvectz(j)
 
        do i=1, nsp
          cij(i)=0.d0
          sik(i)=0.d0
        enddo
 
        do i=1, na
          qtr = qx*fullpos(i,1,k)+qy*fullpos(i,2,k)+qz*fullpos(i,3,k)
          sini = sin(qtr)
          cosi = cos(qtr)
          sp = lot(i)
          sik(sp) = sik(sp) + sini
          cij(sp) = cij(sp) + cosi
        enddo
 
        do i=1, nsp
        do m=1, nsp
          sij(l,i,m) = sij(l,i,m) + cij(i)*cij(m) + sik(i)*sik(m)
        enddo
        enddo
 
     endif
    enddo
  enddo
 
END SUBROUTINE
 
#ifdef OPENMP
!************************************************************
!
! This subroutine computes the sine and cosine sums from the
! configuration for all q-vectors needed
! OpenMP // on Qvect
!
SUBROUTINE fourier_trans_qvect ()
 
  USE parameters
 
  !$ USE OMP_LIB
  IMPLICIT NONE
 
  INTEGER :: numth, sp
  DOUBLE PRECISION :: qx, qy, qz, qtr, sini, cosi
 
  numth = omp_get_max_threads()
  if (number_of_qvect.lt.numth) numth=number_of_qvect
  ! OpemMP on Qvect
  !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
  !$OMP& PRIVATE(qx, qy, qz, cij, sik, qtr, sini, cosi, i, j, k, l, m, sp) &
  !$OMP& SHARED(NUMTH, qvectx, qvecty, qvectz, &
  !$OMP& FULLPOS, NS, NSP, NA, LOT, DELTA_Q, NUMBER_OF_QVECT, NQ, modq, qvmin, Sij)
  !$OMP DO SCHEDULE(STATIC,NUMBER_OF_QVECT/NUMTH)
  do j=1, number_of_qvect
 
    l=anint((modq(j)-qvmin)/delta_q)+1
    if (l .le. nq) then
 
      qx=qvectx(j)
      qy=qvecty(j)
      qz=qvectz(j)
 
      do k=1, ns
 
        do i=1, nsp
          cij(i)=0.d0
          sik(i)=0.d0
        enddo
 
        do i=1, na
          qtr = qx*fullpos(i,1,k)+qy*fullpos(i,2,k)+qz*fullpos(i,3,k)
          sini = sin(qtr)
          cosi = cos(qtr)
          sp = lot(i)
          sik(sp) = sik(sp) + sini
          cij(sp) = cij(sp) + cosi
        enddo
 
        do i=1, nsp
        do m=1, nsp
          !$OMP ATOMIC
          sij(l,i,m) = sij(l,i,m) + cij(i)*cij(m) + sik(i)*sik(m)
        enddo
        enddo
 
      enddo
 
    endif
 
  enddo
 
  !$OMP END DO NOWAIT
  !$OMP END PARALLEL
 
END SUBROUTINE
#endif
 
INTEGER FUNCTION sk_save()
 
USE parameters
 
INTEGER :: nsq
DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: sqtab
 
INTERFACE
  LOGICAL FUNCTION fzbt (NDQ, SQIJ)
    USE parameters
    INTEGER, INTENT(IN) :: ndq
    DOUBLE PRECISION, DIMENSION(NDQ,NSP,NSP), INTENT(IN) :: sqij
  END FUNCTION
END INTERFACE
 
i=0
do j=1, nq
  if (degeneracy(j) .gt. 0) i=i+1
enddo
nsq=i
 
if (nsq .gt. 0) then  ! If wave vectors exist
 
  if (allocated(sqtab)) deallocate(sqtab)
  allocate(sqtab(nsq), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: SK_SAVE"//char(0), "Table: SQTAB"//char(0))
    sk_save = 0
    goto 001
  endif
 
  i = 0;
  do k=1, nq
    if (degeneracy(k).gt.0) then
      i=i+1
      sqtab(i)= k_point(k)
    endif
  enddo
  call save_xsk (nsq, sqtab)
 
  i=0
  do k=1, nq
    if (degeneracy(k).gt.0) then
      i=i+1
      sqtab(i)= s(k)
    endif
  enddo
  call save_curve (nsq, sqtab, 0, idsk)
 
  i=0
  do k=1, nq
    if (degeneracy(k).gt.0) then
      i=i+1
      sqtab(i)= (s(k)-1.0)*k_point(k)
    endif
  enddo
  call save_curve (nsq, sqtab, 2, idsk)
 
  i=0
  do k=1, nq
    if (degeneracy(k).gt.0) then
      i=i+1
      sqtab(i)= xs(k)
    endif
  enddo
  call save_curve (nsq, sqtab, 4, idsk)
 
  i=0
  do k=1, nq
    if (degeneracy(k).gt.0) then
      i=i+1
      sqtab(i)= (xs(k)-1.0)*k_point(k)
    endif
  enddo
  call save_curve (nsq, sqtab, 6, idsk)
 
  l = 8
  do i=1, nsp
    do j=1, nsp
      m=0
      do k=1, nq
        if (degeneracy(k).gt.0) then
          m=m+1
          sqtab(m)=sij(k,i,j)
        endif
      enddo
      call save_curve (nsq, sqtab, l, idsk)
      l=l+2
    enddo
  enddo
 
  !  To compute FZ and BT partials
  if (.not.fzbt(nq, sij)) then
    sk_save = 0
    goto 001
  endif
 
  do i=1, nsp
    do j=1, nsp
      m=0
      do k=1, nq
        if (degeneracy(k).gt.0) then
          m=m+1
          sqtab(m)= fzsij(k,i,j)
        endif
      enddo
      call save_curve (nsq, sqtab, l, idsk)
      l=l+2
    enddo
  enddo
  if (nsp .eq. 2) then
    do i=1, 4
      k=0
      do j=1, nq
        if (degeneracy(j).gt.0) then
          k=k+1
          sqtab(k)= btij(j,i)
        endif
      enddo
      call save_curve (nsq, sqtab, l, idsk)
      l=l+2
    enddo
  endif
 
  sk_save=1
 
endif ! If wave vectors exist
 
001 continue
 
if (allocated(fzsij)) deallocate(fzsij)
if (nsp.eq.2 .and. allocated(btij)) deallocate(btij)
if (allocated(sqtab)) deallocate(sqtab)
 
END FUNCTION
 
END FUNCTION s_of_k
 
INTEGER (KIND=c_int) FUNCTION smooth_and_save (DPOINT, CTS, SFC, IDC, NQPTS, DATS) bind (C,NAME='smooth_and_save_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: idc, nqpts, dats
real(kind=c_double), DIMENSION(NQPTS), INTENT(IN) :: dpoint, cts
real(kind=c_double), INTENT(IN) :: sfc
DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: sqtab
 
INTERFACE
  LOGICAL FUNCTION smooth (TABTOLISS, GTOLISS, DIMTOLISS, SIGMALISS)
    INTEGER, INTENT(IN) :: dimtoliss
    DOUBLE PRECISION, INTENT(IN) :: sigmaliss
    DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMTOLISS) :: gtoliss
    DOUBLE PRECISION, INTENT(INOUT), DIMENSION(DIMTOLISS) :: tabtoliss
  END FUNCTION
END INTERFACE
 
if (allocated(sqtab)) deallocate(sqtab)
 
allocate(sqtab(nqpts), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: smooth_and_save"//char(0), "Table: SQTAB"//char(0))
  smooth_and_save=0
  goto 001
endif
 
do k=1, nqpts
  sqtab(k)=cts(k)
enddo
 
if (.not.smooth(sqtab, dpoint, nqpts, sfc)) then
  smooth_and_save=0
  goto 001
endif
 
call save_curve (nqpts, sqtab, idc, dats)
 
smooth_and_save=1
 
001 continue
 
if (allocated(sqtab)) deallocate(sqtab)
 
END FUNCTION
 
DOUBLE PRECISION FUNCTION fqx(TA, Q)
 
USE mendeleiev
 
INTEGER, INTENT(IN) :: ta
DOUBLE PRECISION, INTENT(IN) :: q
DOUBLE PRECISION :: sinla
DOUBLE PRECISION, PARAMETER :: pi=acos(-1.0)
!
! Acta Cryst. (1968). A24, 321
!
! SINLA = ((sin(THETA)/LAMBDA)**2
! and 2 d Sin(THETA) = LAMBDA
! so 2 d = LAMBDA / Sin(THETA)
! with d = 2 PI / Q we get:
! 4.0 PI / Q =  LAMBDA / Sin(THETA)
sinla=(q/(4.0*pi))**2
 
fqx = a1(ta)*exp(-b1(ta)*sinla) &
    + a2(ta)*exp(-b2(ta)*sinla) &
    + a3(ta)*exp(-b3(ta)*sinla) &
    + a4(ta)*exp(-b4(ta)*sinla) + c(ta)
 
END FUNCTION