spherical.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
!
 
! Bond order parameter from spherical harmonics
 
INTEGER (KIND=c_int) FUNCTION sphericals (MAXL, SPC, GEO, IDC, COOSPH)
 
  USE parameters
#ifdef OPENMP
!$ USE OMP_LIB
#endif
  IMPLICIT NONE
 
  INTEGER (KIND=c_int), INTENT(IN) :: maxl, spc, geo, idc
  INTEGER (KIND=c_int), DIMENSION(NSP), INTENT(IN) :: coosph
  INTEGER :: nspsh
  INTEGER, DIMENSION(:), ALLOCATABLE :: neigh
  LOGICAL :: sphrun
  DOUBLE PRECISION :: xc, yc, zc
  DOUBLE PRECISION :: sr, st, sp
  DOUBLE PRECISION, DIMENSION(0:MAXL,-MAXL:MAXL) :: hsp
! DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: THSP
  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: athsp
  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: sptshp
  DOUBLE PRECISION, DIMENSION(0:MAXL) :: spha
#ifdef OPENMP
  INTEGER :: numth
  LOGICAL :: doatoms
#endif
  INTERFACE
    DOUBLE PRECISION FUNCTION plegendre (l, m, x)
      INTEGER, INTENT(IN) :: l, m
      DOUBLE PRECISION, INTENT(IN) :: x
    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 (allocated(neigh)) deallocate(neigh)
  allocate(neigh(nsp), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: sphericals"//char(0), "Table: NEIGH"//char(0))
    sphericals=0
    goto 001
  endif
 
!  if (allocated(THSP)) deallocate(THSP)
!  allocate(THSP(NBONDS,0:MAXL,-MAXL:MAXL), STAT=ERR)
!  if (ERR .ne. 0) then
!    call show_error ("Impossible to allocate memory"//CHAR(0), &
!                     "Function: sphericals"//CHAR(0), "Table: THSP"//CHAR(0))
!    sphericals=0
!    goto 001
!  endif
  if (allocated(athsp)) deallocate(athsp)
  allocate(athsp(0:maxl,-maxl:maxl), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: sphericals"//char(0), "Table: ATHSP"//char(0))
    sphericals=0
    goto 001
  endif
  if (allocated(sptshp)) deallocate(sptshp)
  allocate(sptshp(0:maxl,-maxl:maxl), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: sphericals"//char(0), "Table: SPTSHP"//char(0))
    sphericals=0
    goto 001
  endif
 
  athsp(:,:)=0.0d0
  sptshp(:,:)=0.0d0
  anbonds=0
  nspsh=0
#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.
 
  if (doatoms) then
    if (na.lt.numth) numth=na
#ifdef DEBUG
    write (6, *) "OpenMP on atoms, NUMTH= ",numth
#endif
    ! OpemMP on atoms
    do i=1, ns
 
      !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
      !$OMP& PRIVATE(SPHRUN, NEIGH, Dij, Rij, j, k, l, m, XC, YC, ZC, SR, ST, SP, HSP) &
      !$OMP& SHARED(NUMTH, i, NA, NSP, NCELLS, LOT, SPC, CONTJ, VOISJ, NSPSH, ATHSP, SPTSHP, ANBONDS, COOSPH, MAXL)
      !$OMP DO SCHEDULE(STATIC,NA/NUMTH)
      do j=1, na
 
        if (lot(j) .eq. spc+1) then
          sphrun=.true.
          neigh(:)=0
          do k=1, contj(j,i)
            neigh(lot(voisj(k,j,i)))=neigh(lot(voisj(k,j,i)))+1
          enddo
 
          !$OMP ATOMIC
          nspsh=nspsh+contj(j,i)
          do k=1, nsp
            if (neigh(k) .ne. coosph(k)) then
              sphrun=.false.
              exit
            endif
          enddo
 
          if (sphrun) then
            !$OMP ATOMIC
            anbonds=anbonds+contj(j,i)
          endif
 
          do k=1, contj(j,i)
            if (ncells .gt. 1) then
              dij = calcdij(rij, j, voisj(k,j,i), i, i, i)
            else
              dij = calcdij(rij, j, voisj(k,j,i), i, i, 1)
            endif
            xc=rij(1)
            yc=rij(2)
            zc=rij(3)
!           NBONDS=NBONDS+1
            call cart2spher(xc, yc, zc, sr, st, sp)
            do l=0, maxl
              do m=0, l
                hsp(l,m) = 0.0d0
                hsp(l,m) = plegendre(l, m, st) * cos(sp * m)
                if (m .ne. 0) hsp(l,-m) = (-1)**m*hsp(l,m)
              enddo
              do m=-l, l
!                THSP(NBONDS,l,m)=HSP(l,m)
                !$OMP ATOMIC
                sptshp(l,m)=sptshp(l,m)+hsp(l,m)
                if (sphrun) then
                  !$OMP ATOMIC
                  athsp(l,m)=athsp(l,m)+hsp(l,m)
                endif
              enddo
            enddo
          enddo
        endif
      enddo
      !$OMP END DO NOWAIT
      !$OMP END PARALLEL
    enddo
 
  else
#ifdef DEBUG
    write (6, *) "OpenMP on MD steps, NUMTH= ",numth
#endif
  ! OpemMP on MD steps
  !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
  !$OMP& PRIVATE(SPHRUN, NEIGH, Dij, Rij, i, j, k, l, m, XC, YC, ZC, SR, ST, SP, HSP) &
  !$OMP& SHARED(NUMTH, NS, NA, NSP, NCELLS, LOT, SPC, CONTJ, VOISJ, NSPSH, ATHSP, SPTSHP, ANBONDS, COOSPH, MAXL)
  !$OMP DO SCHEDULE(STATIC,NS/NUMTH)
#endif
    do i=1, ns
      do j=1, na
 
        if (lot(j) .eq. spc+1) then
          sphrun=.true.
          neigh(:)=0
          do k=1, contj(j,i)
            neigh(lot(voisj(k,j,i)))=neigh(lot(voisj(k,j,i)))+1
          enddo
#ifdef OPENMP
          !$OMP ATOMIC
#endif
          nspsh=nspsh+contj(j,i)
          do k=1, nsp
            if (neigh(k) .ne. coosph(k)) then
              sphrun=.false.
              exit
            endif
          enddo
 
          if (sphrun) then
#ifdef OPENMP
            !$OMP ATOMIC
#endif
            anbonds=anbonds+contj(j,i)
          endif
 
          do k=1, contj(j,i)
            if (ncells .gt. 1) then
              dij = calcdij(rij, j, voisj(k,j,i), i, i, i)
            else
              dij = calcdij(rij, j, voisj(k,j,i), i, i, 1)
            endif
            xc=rij(1)
            yc=rij(2)
            zc=rij(3)
!           NBONDS=NBONDS+1
            call cart2spher(xc, yc, zc, sr, st, sp)
            do l=0, maxl
              do m=0, l
                hsp(l,m) = 0.0d0
                hsp(l,m) = plegendre(l, m, st) * cos(sp * m)
                if (m .ne. 0) hsp(l,-m) = (-1)**m*hsp(l,m)
              enddo
              do m=-l, l
!                THSP(NBONDS,l,m)=HSP(l,m)
#ifdef OPENMP
                !$OMP ATOMIC
#endif
                sptshp(l,m)=sptshp(l,m)+hsp(l,m)
                if (sphrun) then
#ifdef OPENMP
                  !$OMP ATOMIC
#endif
                  athsp(l,m)=athsp(l,m)+hsp(l,m)
                endif
              enddo
            enddo
          enddo
        endif
      enddo
    enddo
#ifdef OPENMP
    !$OMP END DO NOWAIT
    !$OMP END PARALLEL
  endif
#endif
 
  if (geo .eq. 0) then
    spha(:)=0.0d0
    if (nspsh .gt. 0) then
      do l=0, maxl
        do m=-l, l
          spha(l)=spha(l)+(sptshp(l,m)/nspsh)**2
        enddo
        spha(l)=sqrt(4*pi*spha(l)/(2*l+1))
      enddo
    endif
    call save_curve (maxl+1, spha, idc-1, idsp)
  endif
  spha(:)=0.0d0
  if (anbonds .gt. 0) then
    do l=0, maxl
      do m=-l, l
        spha(l)=spha(l)+(athsp(l,m)/anbonds)**2
      enddo
      spha(l)=sqrt(4*pi*spha(l)/(2*l+1))
    enddo
  endif
  call save_curve (maxl+1, spha, idc, idsp)
 
  sphericals = 1
 
  001 continue
 
!  if (allocated(THSP)) deallocate(THSP)
  if (allocated(athsp)) deallocate(athsp)
  if (allocated(sptshp)) deallocate(sptshp)
  if (allocated(neigh)) deallocate(neigh)
 
END FUNCTION
 
SUBROUTINE cart2spher(XC, YC, ZC, RS, TS, PS)
 
DOUBLE PRECISION, INTENT(IN) :: XC, YC, ZC
DOUBLE PRECISION, INTENT(OUT) :: RS, TS, PS
 
rs=0.0d0
ts=0.0d0
ps=0.0d0
 
rs= sqrt(xc**2 + yc**2 + zc**2)
ts= acos(zc/rs)
ps= atan2(xc, yc)
 
END SUBROUTINE
 
! Adapted from Numercial Recipies !
 
DOUBLE PRECISION FUNCTION plegendre (l, m, x)
 
  DOUBLE PRECISION, PARAMETER :: pi=acos(-1.0)
  INTEGER, INTENT(IN) :: l, m
  DOUBLE PRECISION, INTENT(IN) :: x
  DOUBLE PRECISION :: y
  INTEGER :: i, ll
  DOUBLE PRECISION :: fact, oldfact, pll, pmm, pmmp1, omx2
 
  y=cos(x)
  pmm=1.0
  if (m > 0) then
    omx2=(1.0-y)*(1.0+y)
    fact=1.0
    i = 1
    do while (i.le.m)
      pmm = pmm * (omx2*fact/(fact+1.0))
      fact = fact + 2.0
      i = i+1
    enddo
  endif
  pmm=sqrt((2*m+1)*pmm/(4.0*pi))
  if (mod(m,2) .eq. 1) pmm=-pmm
 
  if (l .eq. m) then
    plegendre=pmm
  else
    pmmp1=y*sqrt(2.0*m+3.0)*pmm
    if (l .eq. (m+1)) then
      plegendre=pmmp1
    else
      oldfact=sqrt(2.0*m+3.0)
      ll=m+2
      do while (ll.le.l)
        fact=sqrt((4.0*ll*ll-1.0)/(ll*ll-m*m))
        pll=(y*pmmp1-pmm/oldfact)*fact
        oldfact=fact
        pmm=pmmp1
        pmmp1=pll
        ll=ll+1
      enddo
      plegendre=pll
    endif
  endif
 
END FUNCTION
 
DOUBLE PRECISION FUNCTION plgndr(l, m, x)
 
  DOUBLE PRECISION, PARAMETER :: pi=acos(-1.0)
  INTEGER :: j
  INTEGER, INTENT(IN) :: l, m
  DOUBLE PRECISION, INTENT(IN) :: x
  DOUBLE PRECISION :: prod
 
  INTERFACE
    DOUBLE PRECISION FUNCTION plegendre(l, m, x)
      INTEGER, INTENT(IN) :: l, m
      DOUBLE PRECISION, INTENT(IN) :: x
    END FUNCTION
  END INTERFACE
 
  prod=1.0d0
  j=l-m+1
  do while (j.le.l+m)
    prod = prod*j
    j=j+1
  enddo
 
  plgndr=sqrt(4.0*pi*prod/(2*l+1))*plegendre(l,m,x)
 
END FUNCTION