cqvf.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 cqvf.F90
!! @short q vectors selection for the reciprocal calculation of the S(k)
!! @author Sébastien Le Roux <sebastien.leroux@ipcms.unistra.fr>
 
INTEGER (KIND=c_int) FUNCTION cqvf (QMAX, QMIN, NQ, PROBA, LIMQ) bind (C,NAME='cqvf_')
 
!
! Determines the q-vectors for which the q-dependent
! correlation functions are computed - general cell
!
 
  USE parameters
 
  IMPLICIT NONE
 
  INTEGER (KIND=c_int), INTENT(IN) :: nq
  REAL (kind=c_double), INTENT(IN) :: proba, limq
  REAL (kind=c_double), INTENT(IN) :: qmax, qmin
 
  LOGICAL :: kpts, keepkpts
  INTEGER :: hkl, q_index, seed
  INTEGER :: klow, llow
  INTEGER, DIMENSION(3) :: nkpts
  DOUBLE PRECISION, DIMENSION(3) :: qmrecip
  DOUBLE PRECISION :: qmax2, qmin2, limq2
  DOUBLE PRECISION :: kpx, kpy, kpz, keep
 
  INTERFACE
    DOUBLE PRECISION FUNCTION ran3 (idnum)
      INTEGER, INTENT(IN) :: idnum
    END FUNCTION
  END INTERFACE
 
  qmrecip(:) = 0.0d0
  do i=1, ncells
    do j=1, 3
      qmrecip(j) = qmrecip(j)+the_box(i)%modr(j)
    enddo
  enddo
  qmrecip(:) = qmrecip(:)/ncells
  do i=1, 3
    nkpts(i) = anint(qmax/qmrecip(i))+1
  enddo
  hkl=0
! To simplify the case of the calculation in the case of a simple cubic box
  if (overall_cubic) then
    qmax2=(qmax/qmin)**2
    qmin2=1.0d0
    limq2=(limq/qmin)**2
  else
    qmax2=qmax**2                    ! mod(q_max)**2
    limq2=limq**2
    qmin2=qmin**2                    ! min module of the reciprocal lattice vectors - lattice.f90
  endif
  q_index=0
  qvmax=0.0d0
  qvmin=50.0d0
 
do i=1, 2
 
! The first iteration to evaluate the number of k-points
! to be saved for the analysis, and the second to store them.
 
  if (i .eq. 2) then
 
    number_of_qmod=q_index
    if (allocated(qvectx)) deallocate(qvectx)
    allocate(qvectx(number_of_qmod), stat=err)
    if (err .ne. 0) then
      call show_error ("Impossible to allocate memory"//char(0), &
                       "Function: COMP_Q_VAL_FULL"//char(0), "Table: qvectx"//char(0))
      cqvf=0
      goto 001
    endif
    if (allocated(qvecty)) deallocate(qvecty)
    allocate(qvecty(number_of_qmod), stat=err)
    if (err .ne. 0) then
      call show_error ("Impossible to allocate memory"//char(0), &
                       "Function: COMP_Q_VAL_FULL"//char(0), "Table: qvecty"//char(0))
      cqvf=0
      goto 001
    endif
    if (allocated(qvectz)) deallocate(qvectz)
    allocate(qvectz(number_of_qmod), stat=err)
    if (err .ne. 0) then
      call show_error ("Impossible to allocate memory"//char(0), &
                       "Function: COMP_Q_VAL_FULL"//char(0), "Table: qvectz"//char(0))
      cqvf=0
      goto 001
    endif
    if (allocated(modq)) deallocate(modq)
    allocate(modq(number_of_qmod), stat=err)
    if (err .ne. 0) then
      call show_error ("Impossible to allocate memory"//char(0), &
                       "Function: COMP_Q_VAL_FULL"//char(0), "Table: modq"//char(0))
      cqvf=0
      goto 001
    endif
    q_index=0
 
  endif
 
  do h=0, nkpts(1)
 
    if (h .eq. 0) then
      klow=0
    else
      klow=-nkpts(2)
    endif
 
    do k=klow, nkpts(2)
 
      if (k.eq.0 .and. h.eq.0) then
        llow=0
      else
        llow=-nkpts(3)
      endif
 
      do l=llow, nkpts(3)
 
        if (overall_cubic) then
          kpx= h
          kpy= k
          kpz= l
        else
          ! Note : because we only consider the first step this is an approximation in the case of NPT trajectory
          kpx= 2.0d0*pi*(h*the_box(1)%lrecp(1,1) + k*the_box(1)%lrecp(2,1) + l*the_box(1)%lrecp(3,1))
          kpy= 2.0d0*pi*(h*the_box(1)%lrecp(1,2) + k*the_box(1)%lrecp(2,2) + l*the_box(1)%lrecp(3,2))
          kpz= 2.0d0*pi*(h*the_box(1)%lrecp(1,3) + k*the_box(1)%lrecp(2,3) + l*the_box(1)%lrecp(3,3))
        endif
        qvmod= kpx**2 + kpy**2 + kpz**2
 
!       we want the maximum number of k-points in the FSDP part of the spectra
!       therefore we choose to keep all qvectors with qmod in:
!                   minr < qmod < limq
!       where 'minr' is minimum q modulus accessible considering the analyzed
!       lattice ie. the minimum modulus of the reciprocal cell vectors,
!       'limq' is the limit to be fixed for the FSDP part of the spectra
!       Thereafter: QMIN=minr**2 and LIMQ2=limq**2.
!       The qvectors with a 'qmod > limq' are accepted with a probability of 'PROBA'
 
        if (qvmod .le. qmax2 .and. qvmod .ge. qmin2) then
          kpts=.true.
        else
          kpts=.false.
        endif
        if (i .eq. 1 .and. kpts) then
          q_index=q_index+2
        elseif (i .eq. 2 .and. kpts) then
          if (qvmod .le. limq2) then
            keepkpts=.true.
          else
            seed=173932
            keep = ran3(seed+h**3+k**2+l**5)
            if (keep .le. proba) then
              keepkpts=.true.
            else
              keepkpts=.false.
            endif
          endif
          if (keepkpts) then
            q_index=q_index+1
            qvectx(q_index)=kpx
            qvecty(q_index)=kpy
            qvectz(q_index)=kpz
            modq(q_index)=sqrt(qvmod)
            qvmax=max(qvmax,modq(q_index))
            qvmin=min(qvmin,modq(q_index))
            if (h.ne.0 .or. k.ne.0 .or. l.ne.0) then
              q_index=q_index+1
              qvectx(q_index)=-kpx
              qvecty(q_index)=-kpy
              qvectz(q_index)=-kpz
              modq(q_index)=sqrt(qvmod)
            endif
          endif
        endif
 
      enddo ! l
 
    enddo ! k
 
  enddo ! h
 
enddo
 
number_of_qvect=q_index
 
! NQ is given in input
! the value of each K_POINT is found in the
! interval |qvmax - qvmin| and then we compute
! the degeneracy of each K_POINT in q modulus.
 
if (allocated(degeneracy)) deallocate(degeneracy)
allocate(degeneracy(nq+1), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: COMP_Q_VAL_FULL"//char(0), "Table: degeneracy"//char(0))
  cqvf=0
  goto 001
endif
degeneracy(:)=0
delta_q=(qvmax-qvmin)/nq
 
do i=1, number_of_qvect
  hkl=anint((modq(i)-qvmin)/delta_q)+1
  degeneracy(hkl)=degeneracy(hkl)+1
enddo
 
if (allocated(k_point)) deallocate(k_point)
allocate(k_point(nq), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: COMP_Q_VAL_FULL"//char(0), "Table: K_POINT"//char(0))
  cqvf=0
  goto 001
endif
 
do i=1, nq
  k_point(i)=(i-1.0)*delta_q+qvmin
enddo
 
if (overall_cubic) then
  do i=1, number_of_qvect
    qvectx(i)=qvectx(i)*qmin
    qvecty(i)=qvecty(i)*qmin
    qvectz(i)=qvectz(i)*qmin
  enddo
  do i=1, nq
    k_point(i)=k_point(i)*qmin
  enddo
endif
 
!do i=1, NQ
!  write (6, *) "i= ",i,", K_POINT(i)= ",K_POINT(i)
!enddo
 
cqvf=1
 
001 continue
 
END FUNCTION