gr.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-2024 by CNRS and University of Strasbourg
!
 
INTEGER (KIND=c_int) FUNCTION g_of_r (NDR, DTR, FCR) bind (C,NAME='g_of_r_')
 
! Radial Pair Distribution function
 
!
!             -2                             V
!    g(r)= rho    < Sum Sum d(ri)d(rj-r) > = -- < Sum Sum d(r-rij) >
!                    i  j!i                  N²    i  j!i
!
!
!
!
!                1                    rho(b)                 rho(b)
!    g  (r)= ------------  Sum Sum  ----------  = Sum Sum ------------------
!     ab     4*PI*Delat(r)  a  b!a  d(ab)²*N(a)    a  b!a  Vshell(d(ab))*N(a)
!
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: ndr, fcr
real(kind=c_double), INTENT(IN) :: dtr
DOUBLE PRECISION :: hcap1, hcap2, vcap, dgr
DOUBLE PRECISION :: norm_fact, grlim
DOUBLE PRECISION :: suml, xsuml
LOGICAL :: is_crystal=.false.
#ifdef OPENMP
INTEGER :: numth
LOGICAL :: doatoms
#endif
 
INTERFACE
  LOGICAL FUNCTION allocgr(NDR)
    INTEGER, INTENT(IN) :: ndr
  END FUNCTION
  LOGICAL FUNCTION grbt(GrToBT, NDTR)
    USE parameters
    INTEGER, INTENT(IN) :: ndtr
    DOUBLE PRECISION, DIMENSION(NDTR,NSP,NSP), INTENT(IN) :: grtobt
  END FUNCTION
  DOUBLE PRECISION FUNCTION calcdij (R12, AT1, AT2, STEP_1, STEP_2, SID)
    DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: r12
    INTEGER, INTENT(IN) :: at1, at2, step_1, step_2, sid
  END FUNCTION
END INTERFACE
 
if (.not. allocgr(ndr)) then
  g_of_r = 0
  goto 001
endif
 
do i=1, ndr+1
  shell_vol(i) = 0.0d0
  shell_vol(i) = 4.0d0*pi*(((i-1)*dtr+dtr)**3 - ((i-1)*dtr)**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. i*dtr.gt.mbox) then
    hcap1=i*dtr-mbox
    if ((i-1)*dtr .le. mbox) then
      hcap2=0.0d0
    else
      hcap2=(i-1)*dtr-mbox
    endif
    vcap=hcap1**2*(3*i*dtr - hcap1) - hcap2**2*(3*(i-1)*dtr - hcap2)
    vcap=vcap*pi*2
    shell_vol(i)=shell_vol(i)-vcap
  endif
enddo
 
grlim = ndr*dtr
grlim = grlim*grlim
 
#ifdef OPENMP
  numth = omp_get_max_threads()
  doatoms=.false.
  if (ns.lt.numth) then
    if (numth .ge. 2*(ns-1)) then
      doatoms=.true.
    else
      numth=ns
    endif
  endif
 
  if (all_atoms) doatoms=.true.
#endif
 
if (is_crystal) then
  ! To write the case of highly distoreded crystal
else
#ifdef OPENMP
  if (doatoms) then
    if (na.lt.numth) numth=na
    ! OpemMP on atoms only
    !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
    !$OMP& PRIVATE(Dgr, Dij, Rij, GR_INDEX, NORM_FACT, i, j, k, l, m, n) &
    !$OMP& SHARED(NUMTH, NS, NA, NCELLS, LOT, NBSPBS, SHELL_VOL, DTR, MEANVOL, GRLIM, Gij, Dn, NSP, NDR)
    !$OMP DO SCHEDULE(STATIC,NA-1/NUMTH)
    do i=1, na-1
      do k=1, ns
        do j=i+1, na
          if (ncells .gt. 1) then
            dij = calcdij(rij,i,j,k,k,k)
          else
            dij = calcdij(rij,i,j,k,k,1)
          endif
          if (dij <= grlim) then
            l=lot(i)
            m=lot(j)
            dgr = sqrt(dij)
            if (l .eq. m) then
              n = nbspbs(l)-1
            else
              n = nbspbs(l)
            endif
            gr_index = int(dgr/dtr)+1
            norm_fact = 1.0d0/(shell_vol(gr_index)*dble(n))
            !$OMP ATOMIC
            gij(gr_index,l,m,k) = gij(gr_index,l,m,k) + norm_fact/(dble(nbspbs(m))/meanvol)
            !$OMP ATOMIC
            dn(gr_index,l,m,k) = dn(gr_index,l,m,k) + 1.0d0
          endif
        enddo
      enddo
    enddo
    !$OMP END DO NOWAIT
    !$OMP END PARALLEL
  else
    ! OpemMP on MD steps
    !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
    !$OMP& PRIVATE(Dgr, Dij, Rij, GR_INDEX, NORM_FACT, i, j, k, l, m, n) &
    !$OMP& SHARED(NUMTH, NS, NA, NCELLS, LOT, NBSPBS, SHELL_VOL, DTR, MEANVOL, GRLIM, Gij, Dn, NSP, NDR)
    !$OMP DO SCHEDULE(STATIC,NS/NUMTH)
#endif
    do k=1, ns
      do i=1, na-1
        do j=i+1, na
          if (ncells .gt. 1) then
            dij = calcdij(rij,i,j,k,k,k)
          else
            dij = calcdij(rij,i,j,k,k,1)
          endif
          if (dij <= grlim) then
            l=lot(i)
            m=lot(j)
            dgr = sqrt(dij)
            if (l .eq. m) then
              n = nbspbs(l)-1
            else
              n = nbspbs(l)
            endif
            gr_index = int(dgr/dtr)+1
            norm_fact = 1.0d0/(shell_vol(gr_index)*dble(n))
#ifdef OPENMP
            !$OMP ATOMIC
#endif
            gij(gr_index,l,m,k) = gij(gr_index,l,m,k) + norm_fact/(dble(nbspbs(m))/meanvol)
#ifdef OPENMP
            !$OMP ATOMIC
#endif
            dn(gr_index,l,m,k) = dn(gr_index,l,m,k) + 1.0d0
          endif
        enddo
      enddo
    enddo
#ifdef OPENMP
    !$OMP END DO NOWAIT
    !$OMP END PARALLEL
  endif
#endif
endif
 
do i=1, ndr
  do l=1, ns
    do j=1, nsp
      dn(i,j,j,l) = 2.0d0*dn(i,j,j,l)/dble(nbspbs(j)-1)
    enddo
    do j=1, nsp-1
      do k=j+1, nsp
        dn(i,j,k,l) = dn(i,j,k,l)+dn(i,k,j,l)
        dn(i,k,j,l) = dn(i,j,k,l)
        dn(i,j,k,l) = dn(i,j,k,l)/dble(nbspbs(j))
        dn(i,k,j,l) = dn(i,k,j,l)/dble(nbspbs(k))
      enddo
    enddo
  enddo
enddo
 
do i=2, ndr
  do l=1, ns
    do j=1, nsp
      do k=1, nsp
        dn(i,j,k,l) = dn(i-1,j,k,l) + dn(i,j,k,l)
      enddo
    enddo
  enddo
enddo
 
do i=1, nsp
  do j=1, nsp
    do k=1, ndr
      do l=1, ns
          dn_ij(k,i,j) = dn_ij(k,i,j) + dn(k,i,j,l)
          gr_ij(k,i,j) = gr_ij(k,i,j) + gij(k,i,j,l)
      enddo
      gr_ij(k,i,j) = gr_ij(k,i,j)/dble(ns)
      dn_ij(k,i,j) = dn_ij(k,i,j)/dble(ns)
    enddo
  enddo
enddo
 
do k=1, ndr
  do l=1, nsp
    gr_ij(k,l,l) = 2.0d0*gr_ij(k,l,l)
  enddo
  do j=1, nsp-1
    do i=j+1, nsp
      gr_ij(k,i,j) = gr_ij(k,i,j) + gr_ij(k,j,i)
      gr_ij(k,j,i) = gr_ij(k,i,j)
    enddo
  enddo
enddo
 
if (allocated(dn)) deallocate(dn)
if (allocated(gij)) deallocate(gij)
 
l = 16
do i=1, nsp
  do j=1, nsp
    do k=1, ndr
      grtab(k)=gr_ij(k,i,j)
    enddo
    call save_curve (ndr, grtab, l, idgr)
    l = l+2
    do k=1, ndr
      ggr_ij(k,i,j) = (gr_ij(k,i,j)-1.0)*4.0*pi*(na/meanvol)*(k-0.5)*dtr
      grtab(k) = ggr_ij(k,i,j)
    enddo
    call save_curve (ndr, grtab, l, idgr)
    l = l+2
    do k=1, ndr
      grtab(k)=dn_ij(k,i,j)
    enddo
    call save_curve (ndr, grtab, l, idgr)
    l = l+1
  enddo
enddo
 
if (nsp .eq. 2) then
  if (grbt(gr_ij, ndr)) then
    do k=1, ndr
      grtab(k)=btij(k,1)
    enddo
    call save_curve (ndr, grtab, l, idgr)
    l = l+2
    do k=1, ndr
      grtab(k)=btij(k,2)
    enddo
    call save_curve (ndr, grtab, l, idgr)
    l = l+2
    do k=1, ndr
      grtab(k)=btij(k,3)
    enddo
    call save_curve (ndr, grtab, l, idgr)
  endif
  if (allocated(btij)) deallocate(btij)
endif
 
suml=0.0d0
xsuml=0.0d0
do k=1, nsp
  suml=suml+nscattl(k)*xi(k)
  xsuml=xsuml+xscattl(k)*xi(k)
enddo
suml=suml*suml
xsuml=xsuml*xsuml
 
! Attention bi(NSCATTL(i) dans le code) en fm et bi² en barn = fm*fm*1e-2 = 1e-24cm²
do i=1, ndr
  do j=1, nsp
    do k=1, nsp
      grtot(i) = grtot(i) + xi(j)*xi(k)*nscattl(j)*nscattl(k)*gr_ij(i,j,k)
      ggrtot(i) = ggrtot(i) + xi(j)*xi(k)*nscattl(j)*nscattl(k)*ggr_ij(i,j,k)
      xgrtot(i) = xgrtot(i) + xi(j)*xi(k)*xscattl(j)*xscattl(k)*gr_ij(i,j,k)
      xggrtot(i) = xggrtot(i) + xi(j)*xi(k)*xscattl(j)*xscattl(k)*ggr_ij(i,j,k)
    enddo
  enddo
  grtot(i) = grtot(i)/suml
  drn(i) = ggrtot(i)
  trn(i) =  drn(i) + 4.0*pi*(na/meanvol)*(i-0.5)*dtr*suml
  ggrtot(i) = ggrtot(i)/suml
 
 
  xgrtot(i) = xgrtot(i)/xsuml
  drx(i) = xggrtot(i)
  trx(i) =  drx(i) + 4.0*pi*(na/meanvol)*(i-0.5)*dtr*xsuml
  xggrtot(i) = xggrtot(i)/xsuml
enddo
 
do i=1, ndr
  grtab(i) = grtot(i)
enddo
call save_curve (ndr, grtab, 0, idgr)
 
do i=1, ndr
  grtab(i) = ggrtot(i)
enddo
call save_curve (ndr, grtab, 2, idgr)
 
do i=1, ndr
  grtab(i) = drn(i)
enddo
call save_curve (ndr, grtab, 4, idgr)
 
do i=1, ndr
  grtab(i) = trn(i)
enddo
call save_curve (ndr, grtab, 6, idgr)
 
do i=1, ndr
  grtab(i) = xgrtot(i)
enddo
call save_curve (ndr, grtab, 8, idgr)
 
do i=1, ndr
  grtab(i) = xggrtot(i)
enddo
call save_curve (ndr, grtab, 10, idgr)
 
do i=1, ndr
  grtab(i) = drx(i)
enddo
call save_curve (ndr, grtab, 12, idgr)
 
do i=1, ndr
  grtab(i) = trx(i)
enddo
call save_curve (ndr, grtab, 14, idgr)
 
if (fcr .eq. 1) call fitcutoffs
 
g_of_r=1
 
001 continue
 
if (allocated(dn)) deallocate (dn)
if (allocated(gij)) deallocate (gij)
if (allocated(dn_ij)) deallocate (dn_ij)
if (allocated(grtab)) deallocate (grtab)
if (allocated(ggrtot)) deallocate (ggrtot)
if (allocated(xggrtot)) deallocate (xggrtot)
if (allocated(trn)) deallocate (trn)
if (allocated(trx)) deallocate (trx)
if (allocated(drn)) deallocate (drn)
if (allocated(drx)) deallocate (drx)
if (allocated(ggr_ij)) deallocate (ggr_ij)
if (allocated(gr_ij)) deallocate (gr_ij)
if (allocated(shell_vol)) deallocate(shell_vol)
 
CONTAINS
 
SUBROUTINE fitcutoffs
 
INTERFACE
  INTEGER FUNCTION cutfit (TABTOFIT, NPOINTS)
 
    INTEGER, INTENT(IN) :: npoints
    DOUBLE PRECISION, DIMENSION(NPOINTS), INTENT(IN) :: tabtofit
  END FUNCTION
END INTERFACE
 
if (allocated(gfft)) deallocate(gfft)
allocate(gfft(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Subroutine: FITCUTOFFS"//char(0), "Table: GFFT"//char(0))
endif
 
do l=1, nsp
  do j=1, nsp
    do k=1, ndr
      gfft(k)=gr_ij(k,l,j)
    enddo
    gr_cut(l,j) = (cutfit(gfft, ndr) - 0.5) * dtr
  enddo
enddo
gr_cutoff = (cutfit(grtot, ndr) - 0.5) * dtr
 
call sendcutoffs(nsp, gr_cutoff, gr_cut)
 
if (allocated(gfft)) deallocate(gfft)
 
END SUBROUTINE
 
END FUNCTION
 
INTEGER FUNCTION cutfit (TABTOFIT, NPOINTS)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: npoints
DOUBLE PRECISION, DIMENSION(NPOINTS), INTENT(IN) :: tabtofit
 
o = 1
p = 0
do while (p .eq. 0)
  r=int(o+1+anint(5.0*dble(npoints)/dble(100*ns)))
  if (o.eq.npoints .or. r.gt.npoints) exit
  if (tabtofit(o) .le. tabtofit(o+1)) then
    o=o+1
  else if (tabtofit(o) .le. tabtofit(r)) then
    o=o+1
  else if (tabtofit(o) .le. 0.05) then
    o=o+1
  else
    p=1
  endif
enddo
p = 0
do while (p .eq. 0)
  r=int(o+1+anint(5.0*dble(npoints)/dble(100*ns)))
  if (o.eq.npoints .or. r.gt.npoints) exit
  if (tabtofit(o) .le. tabtofit(o+1)) then
    o=o+1
  else if (tabtofit(o) .le. tabtofit(r)) then
    o=o+1
  else
    p=1
  endif
enddo
 
cutfit = o
 
END FUNCTION
 
LOGICAL FUNCTION allocgr (NDR)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: ndr
 
allocgr=.true.
 
if (allocated(shell_vol)) deallocate(shell_vol)
allocate(shell_vol(ndr+1), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: SHELL_VOL"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(gij)) deallocate(gij)
allocate(gij(ndr+1,nsp,nsp,ns), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Gij"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(dn)) deallocate(dn)
allocate(dn(ndr+1,nsp,nsp,ns), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Dn"//char(0))
  allocgr=.false.
  goto 001
endif
 
gij(:,:,:,:)=0.0d0
dn(:,:,:,:)=0.0d0
 
if (allocated(xgrtot)) deallocate(xgrtot)
allocate(xgrtot(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: XGrTOT"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(grtot)) deallocate(grtot)
allocate(grtot(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: GrTOT"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(gr_ij)) deallocate(gr_ij)
allocate(gr_ij(ndr,nsp,nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Gr_ij"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(dn_ij)) deallocate(dn_ij)
allocate(dn_ij(ndr,nsp,nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Dn_ij"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(grtab)) deallocate(grtab)
allocate(grtab(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: GRTAB"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(ggr_ij)) deallocate(ggr_ij)
allocate(ggr_ij(ndr,nsp,nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Ggr_ij"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(ggrtot)) deallocate(grtot)
allocate(ggrtot(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: GgrTOT"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(xggrtot)) deallocate(xggrtot)
allocate(xggrtot(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: XGgrTOT"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(trn)) deallocate(trn)
allocate(trn(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Trn"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(drn)) deallocate(drn)
allocate(drn(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Drn"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(trx)) deallocate(trx)
allocate(trx(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Trx"//char(0))
  allocgr=.false.
  goto 001
endif
if (allocated(drx)) deallocate(drx)
allocate(drx(ndr), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: ALLOCGR"//char(0), "Table: Drx"//char(0))
  allocgr=.false.
  goto 001
endif
 
grtot(:)=0.0d0
xgrtot(:)=0.0d0
grtab(:)=0.0d0
ggrtot(:)=0.0d0
xggrtot(:)=0.0d0
trn(:)=0.0d0
trx(:)=0.0d0
drn(:)=0.0d0
drx(:)=0.0d0
gr_ij(:,:,:)=0.0d0
dn_ij(:,:,:)=0.0d0
ggr_ij(:,:,:)=0.0d0
 
001 continue
 
END FUNCTION