sq.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 sq.F90
!! @short S(q) analysis: Fourier transform calculation
!! @author Sébastien Le Roux <sebastien.leroux@ipcms.unistra.fr>
 
INTEGER (KIND=c_int) FUNCTION s_of_q (QMAX, QMIN, NQ) bind (C,NAME='s_of_q_')
 
! Total structure factor
! Partials structure factors
 
! Total structure factor:
 
!
!                         00        sin (q*r)
!    S(q)= 1 + 4*PI*rho* /   dr r² ------------ (g(r) - 1)
!                       0              (q*r)
!
 
! or:
 
!             1          Rmax                 sin (q*r)    sin (PI*r/Rmax)
!    S(q)= 1 + 4*PI*rho* /   dr r² (g(r) -1) ----------- * ---------------
!                       0                       (q*r)         PI*r/Rmax
!
 
! Partials structure factors:
 
!
!                                                       Rmax
!    S  (q)=     delta(a,b) + sqrt(x * x ) *4*PI*rho* /   dr r² (g  (r) - 1)
!     ab                            a   b            0            ab
!
!
 
! or:
 
!                                                    Rmax                      sin (q*r)    sin (PI*r/Rmax)
!    S  (q)=  delta(a,b) + sqrt(x * x ) *4*PI*rho*  /   dr r² (g  (r) - go  ) ----------- * ---------------
!     ab                         a   b             0             ab       ab     (q*r)         PI*r/Rmax
!
 
! with x = Nbre(a)/Nbre(tot)
!       a
 
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: nq
real(kind=c_double), INTENT(IN) :: qmax , qmin
DOUBLE PRECISION :: dq
DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: sqtab
 
INTERFACE
  INTEGER FUNCTION recup_data (i, j)
    INTEGER, INTENT(IN) :: i, j
  END FUNCTION
  LOGICAL FUNCTION fzbt (NDQ, SQIJ)
    USE parameters
    INTEGER, INTENT(IN) :: ndq
    DOUBLE PRECISION, DIMENSION(NDQ,NSP,NSP), INTENT(IN) :: sqij
  END FUNCTION
END INTERFACE
 
number_of_qmod = nq
total_density = dble(na)/meanvol
 
if (allocated(q_point)) deallocate(q_point)
allocate(q_point(nq), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_q"//char(0), "Table: Q_POINT"//char(0))
  s_of_q = 0
  goto 001
endif
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_q"//char(0), "Table: S"//char(0))
  s_of_q = 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_q"//char(0), "Table: XS"//char(0))
  s_of_q = 0
  goto 001
endif
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_q"//char(0), "Table: Sij"//char(0))
  s_of_q = 0
  goto 001
endif
 
! 'QMIN' is minimum q modulus accessible considering the analysed
! lattice ie. the minimum modulus of the reciprocal cell vectors,
dq=((qmax-qmin)/dble(nq))
 
s(:)=0.0d0
xs(:)=0.0d0
q_point(:)=0.0d0
do i=1, nq
  q_point(i)= dble(i-1)*dq+qmin
enddo
 
sij(:,:,:)=0.0d0
 
j=0
if (recup_data(j, idgr) .ne. 1) then
  s_of_q = 0
  goto 001
endif
 
j=j+8
if (recup_data(j, idgr) .ne. 1) then
  s_of_q = 0
  goto 001
endif
 
j=j+8
do o=1, nsp
do p=1, nsp
  if (recup_data(j, idgr) .ne. 1) then
    s_of_q = 0
    goto 001
  endif
  j=j+5
enddo
enddo
 
if (.not. fzbt(nq, sij)) then
  s_of_q = 0
  goto 001
endif
 
if (allocated(sqtab)) deallocate(sqtab)
allocate(sqtab(nq), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: s_of_q"//char(0), "Table: SQTAB"//char(0))
  s_of_q = 0
  goto 001
endif
 
l=0
do k=1, nq
  sqtab(k)= s(k)
enddo
call save_curve (nq, sqtab, l, idsq)
l=l+2
 
do k=1, nq
  sqtab(k)= (s(k)-1.0)*q_point(k)
enddo
call save_curve (nq, sqtab, l, idsq)
l=l+2
 
do k=1, nq
  sqtab(k)= xs(k)
enddo
call save_curve (nq, sqtab, l, idsq)
l=l+2
 
do k=1, nq
  sqtab(k)= (xs(k)-1.0)*q_point(k)
enddo
call save_curve (nq, sqtab, l, idsq)
l=l+2
 
do i=1, nsp
do j=1, nsp
  do k=1, nq
    sqtab(k)= sij(k,i,j)
  enddo
  call save_curve (nq, sqtab, l, idsq)
  l=l+2
enddo
enddo
do i=1, nsp
do j=1, nsp
  do k=1, nq
    sqtab(k)= fzsij(k,i,j)
  enddo
  call save_curve (nq, sqtab, l, idsq)
  l=l+2
enddo
enddo
if (nsp .eq. 2) then
  do i=1, 4
    do j=1, nq
      sqtab(j)= btij(j,i)
    enddo
    call save_curve (nq, sqtab, l, idsq)
    l=l+2
  enddo
endif
 
s_of_q = 1
 
001 continue
 
if (allocated(sqtab)) deallocate(sqtab)
if (allocated(q_point)) deallocate(q_point)
if (allocated(sij)) deallocate(sij)
if (allocated(s)) deallocate(s)
if (allocated(xs)) deallocate(xs)
if (allocated(fzsij)) deallocate(fzsij)
if (allocated(btij)) deallocate(btij)
 
END FUNCTION
 
INTEGER (KIND=c_int) FUNCTION send_gr (IC, VAL, DR, RDATA, GDATA) bind (C,NAME='send_gr_')
 
USE parameters
 
INTEGER (KIND=c_int), INTENT(IN) :: ic, val
real(kind=c_double), INTENT(IN) :: dr
real(kind=c_double), DIMENSION(VAL), INTENT(IN) ::  rdata, gdata
DOUBLE PRECISION :: hcap1, hcap2, vcap
INTEGER :: rinit
 
if (allocated(shell_vol)) deallocate(shell_vol)
allocate(shell_vol(val+1), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: send_gr"//char(0), "Table: SHELL_VOL"//char(0))
  send_gr = 0
  goto 001
endif
 
do i=1, val
  shell_vol(i) = 0.0d0
  shell_vol(i) = 4.0d0*pi*((rdata(i)+dr)**3 - (rdata(i)**3))/3
! To take into account the atoms between a/2 and a srqt(3)/2
! We need the small volume they can be found in.
  if (overall_cubic .and. rdata(i)+dr.gt.mbox) then
    hcap1=rdata(i)+dr-mbox
    if (rdata(i) .le. mbox) then
      hcap2=0.0d0
    else
      hcap2=rdata(i)-mbox
    endif
    vcap=hcap1**2*(3*rdata(i) - hcap1) - hcap2**2*(3*rdata(i) - hcap2)
    vcap=vcap*pi*2
    shell_vol(i)=shell_vol(i)-vcap
  endif
enddo
 
j=ic
if (j > 8) then
  m=j-16
  m=m/5
  k=m/nsp+1
  l=m-(k-1)*nsp+1
endif
 
rinit = 1
if (rdata(1) .eq. 0.0) rinit = 2
 
rmax = rdata(val)
do i=1, number_of_qmod
  do n=rinit, val
    phi = q_point(i)*rdata(n)
    fact_rmax = pi*rdata(n)/rmax
!    Sinus_Fact_Rmax = sin(Fact_Rmax)/Fact_Rmax
    sinus_phi = sin(phi)/phi
    if (j .eq. 0) then
      s(i) = s(i) + shell_vol(n)*(gdata(n) - 1)*sinus_phi!*Sinus_Fact_Rmax
    else if (j .eq. 8) then
      xs(i) = xs(i) + shell_vol(n)*(gdata(n) - 1)*sinus_phi
    else
      sij(i,k,l) = sij(i,k,l) + shell_vol(n)*(gdata(n) - 1)*sinus_phi!*Sinus_Fact_Rmax
    endif
  enddo
  if (j .eq. 0) then
    s(i) = 1.0d0 + s(i)*total_density
  else if (j .eq. 8) then
    xs(i) = 1.0d0 + xs(i)*total_density
  else
    if (l .eq. k) then
      sij(i,k,l) = 1.0d0 + xi(l)*sij(i,k,l)*total_density
    else
      sij(i,k,l) = sqrt(xi(k)*xi(l))*sij(i,k,l)*total_density
    endif
  endif
enddo
 
send_gr = 1
 
001 continue
 
if (allocated(shell_vol)) deallocate(shell_vol)
 
END FUNCTION