skt.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 skt.F90
!! @short S(k,t) analysis: dynamic structure factor calculation
!! @author Sébastien Le Roux <sebastien.leroux@ipcms.unistra.fr>
!! @author Noël Jakse <noel.jakse@grenoble-inp.fr>
!
! Notes:
!
! LLM tool (Gemini via Antigravity) was used at few occasions to prepare some sections of this file, including:
!    - To write a first draft version of the 's_of_k_t' function based on my work in the file 'sk.F90'
!    - To write parts of the 'FOURIER_TRANS_QVECT_SKT' routine
!    - To write parts of the 'COMPUTE_SQW' routine
! Overall the inputs provided by the LLM in these occasions were bad,
! But because I used them to get the first raw version of the code,
! I feel compelled to declare it here.
!
INTEGER (KIND=c_int) FUNCTION s_of_k_t (NQ_IN, XA_IN, MIN_IN, N_SETS, SETS_T, &
                                        DELTA_T, Q_NUM, Q_LIST, N_FREQ) bind (C,NAME='s_of_k_t_')
 
! Total and Partial Dynamic Structure Factor Calculation
!
!    S(q,t) = (1/N) * < \sum_{i} \sum_{j} exp( -i * q * ( r_i(t0+t) - r_j(t0) ) ) >
!
! Partial:
!
!    S_ab(q,t) ~ < \rho_a(q, t0+t) * conjg(\rho_b(q, t0)) >
!
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: nq_in  ! Number of delta q
INTEGER (KIND=c_int), INTENT(IN) :: xa_in  ! How to compute X rays
INTEGER (KIND=c_int), INTENT(IN) :: min_in ! Minimum value of correlations
INTEGER (KIND=c_int), INTENT(IN) :: n_sets ! Number of t steps to save, or -1 for all
INTEGER (KIND=c_int), DIMENSION(N_SETS), INTENT(IN) :: sets_t
INTEGER (KIND=c_int), INTENT(IN) :: q_num  ! Number q compute (q,w) data
INTEGER (KIND=c_int), INTENT(IN) :: n_freq ! Number of frequency points
 
real(kind=c_double) :: delta_t
real(kind=c_double), DIMENSION(Q_NUM), INTENT(IN) :: q_list
 
INTEGER :: nsq, pid
INTEGER, DIMENSION(:), ALLOCATABLE :: qid
INTEGER, DIMENSION(:), ALLOCATABLE :: sqw_qlist
DOUBLE PRECISION :: factor, xfactor
DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: sqtab, sqw_tab, sqw_qval
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: skt_tab
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: rho_c, rho_s
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: nsqt, xsqt
DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: localcorr
DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: sqt
 
INTERFACE
  DOUBLE PRECISION FUNCTION fqx(TA, Q)
    INTEGER, INTENT(IN) :: ta
    DOUBLE PRECISION, INTENT(IN) :: q
  END FUNCTION
END INTERFACE
 
allocate(sqt(nq_in, ns-min_in, nsp, nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: SQT"//char(0))
  s_of_k_t = 0
  goto 001
endif
sqt(:,:,:,:) = 0.0d0
 
allocate(nsqt(nq_in, ns-min_in), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: NSQT"//char(0))
  s_of_k_t = 0
  goto 001
endif
nsqt(:,:) = 0.0d0
 
allocate(xsqt(nq_in, ns-min_in), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: XSQT"//char(0))
  s_of_k_t = 0
  goto 001
endif
xsqt(:,:) = 0.0d0
 
! Allocate density arrays RHO_C and RHO_S
ALLOCATE(rho_c(ns, nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: RHO_C"//char(0))
  s_of_k_t = 0
  goto 001
endif
ALLOCATE(rho_s(ns, nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: RHO_S"//char(0))
  s_of_k_t = 0
  goto 001
endif
ALLOCATE(localcorr(ns-min_in, nsp, nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: LocalCorr"//char(0))
  s_of_k_t = 0
  goto 001
endif
 
!t0 = OMP_GET_WTIME ()
call fourier_trans_qvect_skt (min_in) ! Default Q-vector parallelization
!t1 = OMP_GET_WTIME ()
!write (*,*) "temps d’excecution QVT 2:", t1-t0
 
if (allocated(qvectx)) deallocate(qvectx)
if (allocated(qvecty))deallocate(qvecty)
if (allocated(qvectz)) deallocate(qvectz)
if (allocated(modq)) deallocate(modq)
 
! Normalization and weighting (Neutron/X-ray)
 
factor=0.0d0
do i=1, nsp
  factor=factor + nbspbs(i)*nscattl(i)**2
enddo
 
if (xa_in .eq. 1) then
  xfactor=0.0d0
  do i=1, nsp
    xfactor=xfactor + nbspbs(i)*xscattl(i)**2
  enddo
endif
 
do t=1, ns-min_in
 
  do i=1, nq_in
 
    if (degeneracy(i) .gt. 0) then
      do j=1, nsp
        do k=1, nsp
          ! Neutrons
          nsqt(i,t) = nsqt(i,t) + sqt(i,t,j,k) * nscattl(j) * nscattl(k)
          ! X-rays
          if (xa_in .eq. 1) then
            xsqt(i,t) = xsqt(i,t) + sqt(i,t,j,k) * xscattl(j) * xscattl(k)
          else
            ! Use form factors FQX
            xsqt(i,t) = xsqt(i,t) + sqt(i,t,j,k) * fqx(int(xscattl(j)), k_point(i)) * fqx(int(xscattl(k)), k_point(i))
          endif
        enddo
      enddo
 
      ! Normalization
      nsqt(i,t) = nsqt(i,t) / (factor * degeneracy(i))
 
      if (xa_in .eq. 1) then
        xsqt(i,t) = xsqt(i,t) / (xfactor * degeneracy(i))
      else
        ! Function FQX appears in sk.F90
        xfactor = 0.0d0
        do k=1, nsp
          xfactor = xfactor + nbspbs(k) * fqx(int(xscattl(k)), k_point(i))**2
        enddo
        xsqt(i,t) = xsqt(i,t) / (xfactor * degeneracy(i))
      endif
 
      ! Normalize Partials
      do j=1, nsp
        do k=1, nsp
          sqt(i,t,j,k) = sqt(i,t,j,k) / (degeneracy(i) * sqrt(dble(nbspbs(j)*nbspbs(k))))
        enddo
      enddo
 
    endif
 
  enddo
 
enddo
 
if (allocated(sqtab)) deallocate(sqtab)
allocate(sqtab(nq_in), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: SQTAB"//char(0))
  s_of_k_t = 0
  goto 001
endif
 
if (allocated(qid)) deallocate(qid)
allocate(qid(nq_in), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_k_t"//char(0), "Table: QID"//char(0))
  s_of_k_t = 0
  goto 001
endif
sqtab(:)=0.0d0
nsq = 0;
do i=1, nq_in
  if (degeneracy(i) .gt. 0) then
    nsq=nsq+1
    sqtab(nsq)= k_point(i)
    qid(nsq) = i
  endif
enddo
! Saving k-points
call save_xsk (nsq, sqtab)
 
if (skt_save(nsq) .eq. 0) then
  s_of_k_t = 0
  goto 001
endif
 
pid = 8+4*nsp*nsp
if (nsp .eq. 2) pid=pid+8
 
if (n_sets .eq. 1 .and. sets_t(1) .eq. -1) then
  pid = pid*(ns-min_in)
else
  pid = pid*n_sets
endif
 
if (q_num .gt. 0) then
  allocate(sqw_tab(n_freq), stat=err)
  if (err .ne. 0) then
     call show_error ("Impossible to allocate memory"//char(0), &
                      "Function: s_of_k_t"//char(0), "Table: SQW_TAB"//char(0))
     s_of_k_t = 0
     goto 001
  endif
  allocate(sqw_qlist(q_num), stat=err)
  if (err .ne. 0) then
     call show_error ("Impossible to allocate memory"//char(0), &
                      "Function: s_of_k_t"//char(0), "Table: SQW_QLIST"//char(0))
     s_of_k_t = 0
     goto 001
  endif
  allocate(sqw_qval(q_num), stat=err)
  if (err .ne. 0) then
     call show_error ("Impossible to allocate memory"//char(0), &
                      "Function: s_of_k_t"//char(0), "Table: SQW_QVAL"//char(0))
     s_of_k_t = 0
     goto 001
  endif
 
  ! First select all k id for the analysis, as close as possible as the user selection
  ! write (6, *) "Q_NUM= ",Q_NUM
  do i=1, q_num
    j=1
    do while (sqtab(j).lt.q_list(i))
      j=j+1
    enddo
    if (j .gt. 1) then
      if ((sqtab(j) - q_list(i)) .gt. (q_list(i) - sqtab(j-1))) then
        j = j - 1
      endif
    endif
    sqw_qlist(i) = qid(j)
    sqw_qval(i) = k_point(qid(j))
    ! write (6, *) "i= ",i,", Q_LIST(i)= ",Q_LIST(i)," j= ",j,", QID(j)= ",QID(j)," K_POINT(QID(j))= ",K_POINT(QID(j))
  enddo
  call recup_sqw_list (q_num, sqw_qval)
  if (allocated(sqw_qval)) deallocate(sqw_qval)
 
  call compute_sqw (nsqt, q_num, sqw_qlist, pid, 0, 0)
 
  do i=1, nq_in
    do j=1, ns-min_in
      nsqt(i,j) =  (nsqt(i,j)-1.0)*k_point(i)
    enddo
  enddo
  call compute_sqw (nsqt, q_num, sqw_qlist, pid+2, 0, 0)
 
  call compute_sqw (xsqt, q_num, sqw_qlist, pid+4, 0, 0)
 
  do i=1, nq_in
    do j=1, ns-min_in
      xsqt(i,j) =  (xsqt(i,j)-1.0)*k_point(i)
    enddo
  enddo
  call compute_sqw (xsqt, q_num, sqw_qlist, pid+6, 0, 0)
 
  allocate(skt_tab(nq_in, ns-min_in), stat=err)
  if (err .ne. 0) then
     call show_error ("Impossible to allocate memory"//char(0), &
                      "Function: s_of_k_t"//char(0), "Table: SKT_TAB"//char(0))
     s_of_k_t = 0
     goto 001
  endif
 
  pid = pid+8
  do j=1, nsp
    do k=1, nsp
      do l=1, ns-min_in
        do m=1, nq_in
          skt_tab(m,l) = sqt(m,l,j,k)
        enddo
      enddo
      call compute_sqw (skt_tab, q_num, sqw_qlist, pid, j, k)
      pid = pid + 2
    enddo
  enddo
 
endif
 
s_of_k_t = 1
 
001 continue
 
if (allocated(sqtab)) deallocate(sqtab)
if (allocated(skt_tab)) deallocate(skt_tab)
if (allocated(sqw_tab)) deallocate(sqw_tab)
if (allocated(sqw_qlist)) deallocate(sqw_qlist)
if (allocated(sqt)) deallocate(sqt)
if (allocated(nsqt)) deallocate(nsqt)
if (allocated(xsqt)) deallocate(xsqt)
if (allocated(rho_c)) deallocate(rho_c)
if (allocated(rho_s)) deallocate(rho_s)
if (allocated(localcorr)) deallocate(localcorr)
 
CONTAINS
 
!************************************************************
!
! Compute S(q,t) loops over Q-vectors
! OpenMP // on Qvect
!
SUBROUTINE fourier_trans_qvect_skt (MIN_IN)
 
  USE parameters
 
  IMPLICIT NONE
 
  INTEGER, INTENT(IN) :: min_in
 
  INTEGER :: q, numcorr, t_n
  DOUBLE PRECISION :: qx, qy, qz, qtr
  DOUBLE PRECISION :: corr
 
#ifdef OPENMP
  INTEGER :: numth
  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, qtr, i, j, k, l, m, n, q, t) &
  !$OMP& PRIVATE (t_n, NumCorr, RHO_C, RHO_S, LocalCorr, Corr) &
  !$OMP& SHARED(NUMTH, NUMBER_OF_QVECT, SQT, NQ_IN, modq, qvmin, DELTA_Q) &
  !$OMP& SHARED(qvectx, qvecty, qvectz, FULLPOS, NS, NSP, NA, LOT, MIN_IN)
  !$OMP DO SCHEDULE(STATIC,NUMBER_OF_QVECT/NUMTH)
#endif
  do q=1, number_of_qvect
 
    l=anint((modq(q)-qvmin)/delta_q)+1
    if (l .le. nq_in) then
 
      rho_c(:,:) = 0.0d0
      rho_s(:,:) = 0.0d0
      localcorr(:,:,:) = 0.0d0
 
      qx=qvectx(q)
      qy=qvecty(q)
      qz=qvectz(q)
 
      do t=1, ns
        ! Compute density history for this Q vector
        do i=1, na
          j = lot(i)
          qtr = qx*fullpos(i,1,t) + qy*fullpos(i,2,t) + qz*fullpos(i,3,t)
          rho_c(t, j) = rho_c(t, j) + cos(qtr)
          rho_s(t, j) = rho_s(t, j) + sin(qtr)
        enddo
      enddo
 
      ! If 'MIN_IN = 0' and 't = 0', then S(t=0) is the static structure factor
      do t=0, ns-1-min_in
        numcorr = ns-t-min_in
        do t_n=1, numcorr
          do m=1, nsp
            do n=1, nsp
              corr = rho_c(t+t_n, m) * rho_c(t_n, n) + rho_s(t+t_n, m) * rho_s(t_n, n)
              localcorr(t+1, m, n) = localcorr(t+1, m, n) + corr
            enddo
          enddo
        enddo
        ! Normalize by NumCorr here
        localcorr(t+1, :, :) = localcorr(t+1, :, :) / dble(numcorr)
      enddo
 
#ifdef OPENMP
     !$OMP CRITICAL
#endif
      sqt(l, :, :, :) = sqt(l, :, :, :) + localcorr(:, :, :)
#ifdef OPENMP
     !$OMP END CRITICAL
#endif
    endif
 
  enddo
#ifdef OPENMP
  !$OMP END DO
  !$OMP END PARALLEL
#endif
 
END SUBROUTINE
 
!************************************************************
!
! Compute S(q,w) loops over frequencies
!
SUBROUTINE compute_sqw (SKT_TAB, Q_NUM, Q_LIST, PIC, SPA, SPB)
 
  USE parameters
 
  IMPLICIT NONE
 
  INTEGER, INTENT(IN) :: q_num, pic, spa, spb
  INTEGER, DIMENSION(Q_NUM), INTENT(IN) :: q_list
  DOUBLE PRECISION, DIMENSION (NQ_IN,NS-MIN_IN), INTENT(IN) :: skt_tab
 
  INTEGER :: qid, id_q_num, freq
  INTEGER :: shift, cid
  DOUBLE PRECISION :: omega, max_omega, delta_omega, time_val, sqw_val
 
  max_omega = pi / delta_t
  ! delta_omega = max_omega / DBLE(N_FREQ)
  delta_omega = max_omega / dble(n_freq-1)
 
  shift = 8+4*nsp*nsp
  if (nsp .eq. 2) shift=shift+8
  cid = pic
 
  do qid = 1, q_num ! For all selected q points
 
    id_q_num = q_list(qid) ! Select the q point ID number as referenced previously
    do freq = 1, n_freq
 
      omega = dble(freq-1) * delta_omega
      sqw_val = 0.5d0 * skt_tab(id_q_num, 1)
 
      do t = 2, ns-min_in-1
        time_val = dble(t-1) * delta_t
        sqw_val = sqw_val + skt_tab(id_q_num, t) * cos(omega * time_val)
      enddo
 
      time_val = dble(ns-min_in-1) * delta_t
      sqw_val = sqw_val + 0.5d0 * skt_tab(id_q_num, ns-min_in) * cos(omega * time_val)
 
      sqw_val = 2.0d0 * sqw_val * delta_t / pi
      sqw_tab(freq) = sqw_val
 
    enddo
 
    ! Save SQW_TAB here !
    call save_curve (n_freq, sqw_tab, cid, idskt)
 
    if (spa .gt. 0 .and. spb .gt. 0) then
 
      ! For partials only evaluates Faber-Ziman formalism
      sqw_val = 0.0d0
      do freq = 1, n_freq
 
        omega = dble(freq-1) * delta_omega
        if (spa .eq. spb) then
          sqw_val = 0.5d0 * (1.0d0 + (skt_tab(id_q_num, 1) - 1.0d0)/xi(spa))
        else
          sqw_val = 0.5d0 * (1.0d0 + skt_tab(id_q_num, 1)/sqrt(xi(spa)*xi(spb)))
        endif
 
        do t = 2, ns-min_in-1
          time_val = dble(t-1) * delta_t
          if (spa .eq. spb) then
            sqw_val = sqw_val + (1.0d0 + (skt_tab(id_q_num, t) - 1.0d0)/xi(spa)) * cos(omega * time_val)
          else
            sqw_val = sqw_val + (1.0d0 + skt_tab(id_q_num, t)/sqrt(xi(spa)*xi(spb))) * cos(omega * time_val)
          endif
        enddo
 
        time_val = dble(ns-min_in-1) * delta_t
        if (spa .eq. spb) then
          sqw_val = sqw_val + 0.5d0 * (1.0d0 + (skt_tab(id_q_num, t) - 1.0d0)/xi(spa)) * cos(omega * time_val)
        else
          sqw_val = sqw_val + 0.5d0 * (1.0d0 + skt_tab(id_q_num, t)/sqrt(xi(spa)*xi(spb))) * cos(omega * time_val)
        endif
        sqw_val = 2.0d0 * sqw_val * delta_t  / pi
        sqw_tab(freq) = sqw_val
 
      enddo
 
      ! Save SQW_TAB here !
      call save_curve (n_freq, sqw_tab, cid+2*nsp*nsp, idskt)
 
    endif
 
    cid = cid + shift
 
  enddo
 
 
END SUBROUTINE
 
INTEGER FUNCTION skt_save (NSQ)
 
USE parameters
USE mendeleiev
 
INTEGER, INTENT(IN) :: nsq
 
INTEGER :: ndt, tps, cid
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
 
h = 8+4*nsp*nsp
if (nsp .eq. 2) h=h+8
 
if (allocated(sqtab)) deallocate(sqtab)
allocate(sqtab(nq_in), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: SKT_SAVE"//char(0), "Table: SQTAB"//char(0))
  skt_save = 0
  goto 001
endif
if(allocated(sij)) deallocate(sij)
allocate(sij(nq_in,nsp,nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: SKT_SAVE"//char(0), "Table: Sij"//char(0))
  skt_save = 0
  goto 001
endif
 
if (nsq .gt. 0) then  ! If wave vectors exist
 
  sqtab(:)=0.0d0
 
  if (n_sets .eq. 1 .and. sets_t(1) .eq. -1) then
    ndt = ns-min_in
  else
    ndt = n_sets
  endif
 
  do t=1, ndt
 
    cid = (t-1)*sknum
    if (ndt .eq. ns-min_in) then
      tps = t
    else
      tps = sets_t(t)
    endif
 
    i=0
    do k=1, nq_in
      if (degeneracy(k) .gt. 0) then
        i=i+1
        sqtab(i)= nsqt(k,tps)
      endif
    enddo
    call save_curve (nsq, sqtab, cid, idskt)
 
    i=0
    do k=1, nq_in
      if (degeneracy(k) .gt. 0) then
        i=i+1
        sqtab(i)= (nsqt(k,tps)-1.0)*k_point(k)
      endif
    enddo
    call save_curve (nsq, sqtab, cid + 2, idskt)
 
    i=0
    do k=1, nq_in
      if (degeneracy(k) .gt. 0) then
        i=i+1
        sqtab(i)= xsqt(k,tps)
      endif
    enddo
    call save_curve (nsq, sqtab, cid + 4, idskt)
 
    i=0
    do k=1, nq_in
      if (degeneracy(k) .gt. 0) then
        i=i+1
        sqtab(i)= (xsqt(k,tps)-1.0)*k_point(k)
      endif
    enddo
    call save_curve (nsq, sqtab, cid + 6, idskt)
 
    sqtab(:)=0.0d0
    sij(:,:,:)=0.0d0
    l = 8
    do i=1, nsp
      do j=1, nsp
        m=0
 
        do k=1, nq_in
          sij(k,i,j) = sqt(k,tps,i,j)
          if (degeneracy(k) .gt. 0) then
            m=m+1
            sqtab(m)=sij(k,i,j)
          endif
        enddo
        call save_curve (nsq, sqtab, cid + l, idskt)
        l=l+2
      enddo
    enddo
 
  !  To compute FZ and BT partials
    if (.not.fzbt(nq_in, sij)) then
      skt_save = 0
      goto 001
    endif
 
    do i=1, nsp
      do j=1, nsp
        m=0
        do k=1, nq_in
          if (degeneracy(k) .gt. 0) then
            m=m+1
            sqtab(m)= fzsij(k,i,j)
          endif
        enddo
        call save_curve (nsq, sqtab, cid + l, idskt)
        l=l+2
      enddo
    enddo
    if (nsp .eq. 2) then
      do i=1, 4
        k=0
        do j=1, nq_in
          if (degeneracy(j) .gt. 0) then
            k=k+1
            sqtab(k)= btij(j,i)
          endif
        enddo
        call save_curve (nsq, sqtab, cid + l, idskt)
        l=l+2
      enddo
    endif
 
  enddo
 
endif ! If wave vectors exist
 
skt_save=1
 
001 continue
 
if (allocated(fzsij)) deallocate(fzsij)
if (nsp.eq.2 .and. allocated(btij)) deallocate(btij)
if (allocated(sqtab)) deallocate(sqtab)
if(allocated(sij)) deallocate(sij)
if(allocated(fzsij)) deallocate(fzsij)
if(allocated(btij)) deallocate(btij)
 
END FUNCTION
 
END FUNCTION s_of_k_t