bonds.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 bonding (scf, sbf, adv, bdist, bmin, delt_ij, sfil) bind (C,NAME='bonding_')
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: scf, sbf, adv, bdist
real(kind=c_double), INTENT(IN) :: bmin, delt_ij
CHARACTER (KIND=c_char), DIMENSION(*), INTENT(IN) :: sfil
INTEGER, DIMENSION(:), ALLOCATABLE :: gesp
CHARACTER (LEN=scf) :: sfile
DOUBLE PRECISION :: dbd
DOUBLE PRECISION, DIMENSION(3) :: rbd
#ifdef OPENMP
INTEGER :: numth
LOGICAL :: doatoms
#endif
TYPE geometry
  INTEGER :: indice
  INTEGER :: coord
  INTEGER :: ncoord
  INTEGER, DIMENSION(:), ALLOCATABLE :: geo
  TYPE (geometry), POINTER :: first
  TYPE (geometry), POINTER :: last
  TYPE (geometry), POINTER :: next
  TYPE (geometry), POINTER :: prev
END TYPE geometry
 
TYPE(geometry), DIMENSION(:), POINTER :: geot, geop
TYPE(geometry), POINTER :: ga
 
INTERFACE
  LOGICAL FUNCTION allocbonds (alloc)
    !INTEGER, INTENT(IN) :: adv
    LOGICAL, INTENT(IN) :: alloc
  END FUNCTION
    LOGICAL FUNCTION allocedco (alloc)
    LOGICAL, INTENT(IN) :: alloc
  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
  LOGICAL FUNCTION eescs ()
  END FUNCTION
  INTEGER FUNCTION molecules()
  END FUNCTION
END INTERFACE
 
if (adv .eq. 1) then
  if (allocated(statbd)) deallocate(statbd)
  allocate(statbd(nsp,nsp,0:bdist), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: bonding"//char(0), "Table: STATBD"//char(0))
    bonding=0
    goto 001
  endif
  statbd(:,:,:)=0
endif
 
if (.not. allocbonds(.true.)) then
  bonding=0
  goto 001
endif
 
if (allocated(gesp)) deallocate(gesp)
allocate(gesp(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bonding"//char(0), "Table: GESP"//char(0))
  bonding=0
  goto 001
endif
if (allocated(lt_geom)) deallocate(lt_geom)
allocate(lt_geom(nsp,na), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bonding"//char(0), "Table: LT_GEOM"//char(0))
  bonding=0
  goto 001
endif
lt_geom(:,:) = 0
allocate(geot(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bonding"//char(0), "Pointer: GEOT"//char(0))
  bonding=0
  goto 001
endif
allocate(geop(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bonding"//char(0), "Pointer: GEOP"//char(0))
  bonding=0
  goto 001
endif
allocate(ga, stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bonding"//char(0), "Pointer: GA"//char(0))
  bonding=0
  goto 001
endif
do i=1, nsp
  geot(i)%INDICE = 0
  nullify(geot(i)%NEXT)
  nullify(geot(i)%PREV)
  nullify(geot(i)%FIRST)
  nullify(geot(i)%LAST)
  geop(i)%INDICE = 0
  nullify(geop(i)%NEXT)
  nullify(geop(i)%PREV)
  nullify(geop(i)%FIRST)
  nullify(geop(i)%LAST)
enddo
 
! Décompte des nombres de coordination
! Evaluation of the coordination numbers
#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
  do i=1, ns
    ! OpemMP on Atoms
    !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
    !$OMP& PRIVATE(j, l, k, m, n, o, p, DBD, RBD, GESP, GA) &
    !$OMP& SHARED(NUMTH, NS, i, NA, NCELLS, LOT, CONTJ, VOISJ, LA_COUNT, STATBD, adv, bmin, delt_ij)
    !$OMP DO SCHEDULE(STATIC,NA/NUMTH)
    do j=1, na
      k = lot(j)
      l = contj(j,i)
      do m=1, l
        n = voisj(m,j,i)
        o = lot(n)
        la_count(j,o,i)=la_count(j,o,i)+1
        if (adv .eq. 1) then
          if (ncells .gt. 1) then
            dbd = calcdij(rbd, j, n, i, i, i)
          else
            dbd = calcdij(rbd, j, n, i, i, 1)
          endif
          dbd = sqrt(dbd)
          p =int((dbd-bmin)/delt_ij)
          !$OMP ATOMIC
          statbd(k,o,p)=statbd(k,o,p) + 1
        endif
      enddo
    enddo
    !$OMP END DO NOWAIT
    !$OMP END PARALLEL
  enddo
else
  ! OpemMP on MD steps
  !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
  !$OMP& PRIVATE(i, j, k, l, m, n, o, p, DBD, RBD, GESP, GA) &
  !$OMP& SHARED(NUMTH, NS, NA, NCELLS, LOT, CONTJ, VOISJ, LA_COUNT, STATBD, adv, bmin, delt_ij)
  !$OMP DO SCHEDULE(STATIC,NS/NUMTH)
#endif
 do i=1, ns
    do j=1, na
      k = lot(j)
      l = contj(j,i)
      do m=1, l
        n = voisj(m,j,i)
        o = lot(n)
        la_count(j,o,i)=la_count(j,o,i)+1
        if (adv .eq. 1) then
          if (ncells .gt. 1) then
            dbd = calcdij(rbd, j, n, i, i, i)
          else
            dbd = calcdij(rbd, j, n, i, i, 1)
          endif
          dbd = sqrt(dbd)
          p =int((dbd-bmin)/delt_ij)
          !$OMP ATOMIC
          statbd(k,o,p)=statbd(k,o,p) + 1
        endif
      enddo
    enddo
  enddo
#ifdef OPENMP
  !$OMP END DO NOWAIT
  !$OMP END PARALLEL
endif
#endif
 
if (sbf .eq. 1 .and. scf.gt.0) then
  do i=1, scf
    sfile(i:i) = sfil(i)
  enddo
  open(unit=100, file=sfile, action="write", status='unknown')
  do i=1, ns
    write (100, °'("Configuration N",i6)') i
    write (100, *)
    do j=1, na
      write (100, '(A2,i6)') tl(lot(j)), j
      write (100, '(4x,"Nc[tot]= ",i2)') contj(j,i)
      do k=1, nsp
        write (100, '(5x,"Nc[",A2,"]= ",i2)') tl(k), la_count(j,k,i)
        do l=1, contj(j,i)
          if (lot(voisj(l,j,i)) .eq. l) then
            if (ncells .gt. 1) then
              dbd = calcdij(rbd, j, voisj(l,j,i), i, i, i)
            else
              dbd = calcdij(rbd, j, voisj(l,j,i), i, i, 1)
            endif
            write (100, Å'(10x,A2,1x,i6,1x,"at",1x,f7.5,1x,"")') tl(k), voisj(l,j,i), sqrt(dbd)
          endif
        enddo
      enddo
    enddo
  enddo
  close(100)
endif
 
do i=1, ns
  do j=1, na
    k = lot(j)
    l = 0
    do m=1, nsp
      n = la_count(j,m,i)
      l = l + n
      gesp(m) = n
      ma_count(k,m)=ma_count(k,m)+n
      sa_count(k)=sa_count(k)+n
    enddo
    ga => geop(k)
    m = newcoord(1, l, ga, nsp, gesp)
    ga => geot(k)
    n = newcoord(0, l, ga, nsp, gesp)
    if (m > 0 .and. n > 0) then
      lt_geom(k,n) = 1
      togl((i-1)*na+j) = m
      tigl((i-1)*na+j) = n
    else
      bonding=0
      goto 001
    endif
  enddo
enddo
 
n=0
m=0
 
do i=1, nsp
  j = geop(i)%LAST%INDICE
  n = n + j
  ga => geop(i)%FIRST
  do k=1, j
    m = max(m, ga%COORD)
    if (k < j) ga => ga%NEXT
  enddo
enddo
 
call send_coord_opengl (1, ns*na, 0, 0, n, togl)
if (allocated(togl)) deallocate(togl)
 
o = 0
n = 0
m = 20
do i=1, nsp
  j = geot(i)%LAST%INDICE
  o = o + j
  ga => geot(i)%FIRST
  do k=1, j
    n = max(n, ga%COORD)
    m = min(m, ga%COORD)
    if (k < j) ga => ga%NEXT
  enddo
enddo
 
call send_coord_opengl (0, ns*na, m, n, o, tigl)
 
do i=1, nsp
  j = geot(i)%LAST%INDICE
  if (allocated(lgsa)) deallocate(lgsa)
  allocate(lgsa(j), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: bonding"//char(0), "Table: LGSA"//char(0))
    bonding=0
    goto 001
  endif
  ga => geot(i)%FIRST
  do k=1, j
    lgsa(k) = ga%COORD
    if (k < j) ga => ga%NEXT
  enddo
  call init_menu_coordinations (0, i-1, j, lgsa)
enddo
 
if (allocated(lgsa)) deallocate(lgsa)
if (allocated(tigl)) deallocate(tigl)
do i=1, nsp
  sa_count(i)=sa_count(i)/ns/nbspbs(i)
  do j=1, nsp
    ma_count(i,j)=ma_count(i,j)/ns/nbspbs(i)
  enddo
enddo
 
do i=1, nsp
  do j=1, nsp
    mac(j)=ma_count(i,j)
  enddo
  call coordout (i-1, sa_count(i), mac, geop(i)%LAST%INDICE)
enddo
 
do i=1, nsp
  j = geop(i)%LAST%INDICE
  if (allocated(ngsa)) deallocate(ngsa)
  allocate(ngsa(j), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: bonding"//char(0), "Table: NGSA"//char(0))
    bonding=0
    goto 001
  endif
  if (allocated(lgsa)) deallocate(lgsa)
  allocate(lgsa(j), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: bonding"//char(0), "Table: LGSA"//char(0))
    bonding=0
    goto 001
  endif
  ga => geop(i)%FIRST
  call allocate_partial_geo (i-1, j)
  do k=1, j
    lgsa(k)=ga%COORD
    ngsa(k)=ga%NCOORD
    do l=1, nsp
      gesp(l) = ga%GEO(l)
    enddo
    call partial_geo_out (i-1, k-1, nsp, gesp)
    if (k .lt. j) ga => ga%NEXT
  enddo
  call envout (i-1, j, ngsa)
  call init_menu_coordinations (1, i-1, j, lgsa)
  if (allocated(ngsa)) deallocate(ngsa)
  if (allocated(lgsa)) deallocate(lgsa)
enddo
call clean_geom ()
if (allocated(gesp)) deallocate(gesp)
 
if (adv .eq. 1) then
  if (.not.eescs()) then
    bonding=0
    goto 001
  endif
  do i=1, nsp
    do j=1, nsp
      eabl(j)=edgeta(i,j)
      cabl(j)=corta(i,j)
      dabl(j)=defta(i,j)
      etabl(j)=etypea(i,j)
      ctabl(j)=ctypea(i,j)
      dtabl(j)=detypea(i,j)
      tdtabl(j)=tda(i,j)
      ectabl(j)=etda(i,j)
    enddo
    call tetraout (i-1, eabl, cabl, dabl, etabl, ctabl, dtabl, tdtabl, ectabl)
  enddo
endif
if (.not.allocedco(.false.)) then
  bonding=0
  goto 001
endif
 
if (adv .eq. 1) then
 
  if (allocated(eabl)) deallocate(eabl)
  allocate(eabl(0:bdist), stat=err)
  if (err .ne. 0) then
    call show_error ("Impossible to allocate memory"//char(0), &
                     "Function: bonding"//char(0), "Table: EABL"//char(0))
    bonding=0
    goto 001
  endif
  do k=0, bdist
    eabl(k) = 0.0
  enddo
  l=0
  do i=1, nsp
  do j=1, nsp
    m=0
    do k=0, bdist
      m = m + statbd(i,j,k)
    enddo
    if (m .ne. 0) then
      do k=0, bdist
        eabl(k) = 100.0*statbd(i,j,k)/m
      enddo
      call save_curve (bdist+1, eabl, l, idbd)
    else
      call save_curve (0, eabl, l, idbd)
    endif
    l=l+1
  enddo
  enddo
  if (allocated(statbd)) deallocate(statbd)
endif
 
if (.not.allocbonds(.false.)) then
  bonding=0
  goto 001
endif
 
bonding=1
 
001 continue
 
CONTAINS
 
INTEGER FUNCTION newcoord (ID, GEO, GP, NP, GSP)
 
INTEGER, INTENT(IN) :: id, geo, np
INTEGER, DIMENSION(NP), INTENT(IN) :: gsp
TYPE(geometry), INTENT(INOUT), POINTER :: gp
INTEGER :: ab, ac, ad
LOGICAL :: newgeo=.false.
TYPE(geometry), POINTER :: gm
 
if (gp%INDICE .eq. 0) then
  gp%INDICE = 1
  gp%COORD = geo
  gp%NCOORD = 1
  gp%LAST => gp
  gp%FIRST => gp
  newcoord=gp%INDICE
  if (id .eq. 1) then
    allocate(gp%GEO(np))
    do ab=1, np
      gp%GEO(ab) = gsp(ab)
    enddo
  endif
else
 
  gp => gp%FIRST
  ac = gp%LAST%INDICE
  newgeo=.true.
  if (id .eq. 0) then
    do ab=1, ac
      if (gp%COORD .eq. geo) then
        newgeo=.false.
        newcoord = gp%INDICE
        exit
      endif
      if (ab .lt. ac) gp => gp%NEXT
    enddo
  else
    do ab=1, ac
      if (gp%COORD .eq. geo) then
        newgeo=.false.
        do ad=1, np
          if (gp%GEO(ad) .ne. gsp(ad)) then
            newgeo=.true.
            exit
          endif
        enddo
        if (.not.newgeo) exit
      endif
      if (ab .lt. ac) gp => gp%NEXT
    enddo
  endif
  if (newgeo) then
    allocate(gm, stat=err)
    if (err .ne. 0) then
      call show_error ("Impossible to allocate memory"//char(0), &
                       "Function: NEWCOORD"//char(0), "Pointer: GM"//char(0))
      newcoord=0
      goto 001
    endif
    nullify(gm%FIRST)
    nullify(gm%NEXT)
    nullify(gm%LAST)
    nullify(gm%PREV)
    gm%FIRST => gp%FIRST
    gp%NEXT => gm
    gp%FIRST%LAST => gm
    gp%LAST => gm
    gm%LAST => gm
    gm%PREV => gp
    gm%INDICE = gp%INDICE+1
    gm%COORD = geo
    gm%NCOORD = 1
    newcoord = gm%INDICE
    if (id .eq. 1) then
      allocate(gm%GEO(np))
      do ab=1, np
        gm%GEO(ab) = gsp(ab)
      enddo
    endif
    gp => gm
  else
    newcoord = gp%INDICE
    gp%NCOORD = gp%NCOORD+1
  endif
 
endif
 
001 continue
 
END FUNCTION
 
SUBROUTINE clean_geom ()
 
INTEGER :: ab, ac, ad
 
do ab=1, nsp
  ga => geot(ab)%LAST
  ac=ga%INDICE
  do ad=1, ac-1
    ga => ga%PREV
    deallocate (ga%NEXT)
  enddo
enddo
deallocate(geot)
 
do ab=1, nsp
  ga => geop(ab)%LAST
  ac=ga%INDICE
  do ad=1, ac-1
    ga => ga%PREV
    deallocate (ga%NEXT)
  enddo
enddo
deallocate(geop)
 
END SUBROUTINE
 
END FUNCTION
 
SUBROUTINE sendcuts (cspa, cspb, cutab) bind (C,NAME='sendcuts_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: cspa, cspb
real(kind=c_double), INTENT(IN) :: cutab
 
if (cspa .eq. nsp) then
  gr_cutoff=cutab**2
else
  gr_cut(cspa+1,cspb+1) = cutab**2
endif
 
END SUBROUTINE
 
REAL (KIND=c_double) FUNCTION fdmax(use_pbc) bind (C,NAME='fdmax_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: use_pbc
 
if (use_pbc .eq. 0) then
  do l=1,3
    pmin(l) = fullpos(1,l,1)
    pmax(l) = pmin(l)
  enddo
  do i=1, ns
    if (i .eq. 1) then
      k = 2
    else
      k = 1
    endif
    do j=k, na
      do l=1,3
        pmin(l) = min(pmin(l), fullpos(j,l,i))
        pmax(l) = max(pmax(l), fullpos(j,l,i))
      enddo
    enddo
  enddo
  fdmax = 0.0d0
  do l=1,3
    fdmax = fdmax + (pmax(l)-pmin(l))**2
  enddo
  fdmax = sqrt(fdmax)
else
  if (the_box(1)%CUBIC) then
    fdmax = sqrt(2.0)*the_box(1)%minv/2.0
  else
    fdmax = the_box(1)%minv/2.0
  endif
  if (ncells .gt. 1) then
    do l=2, ncells
      if (the_box(l)%CUBIC) then
        fdmax = min(fdmax,sqrt(2.0)*the_box(l)%minv/2.0)
      else
        fdmax = min(fdmax, the_box(l)%minv/2.0)
      endif
    enddo
  endif
 
endif
 
END FUNCTION
 
REAL (KIND=c_double) FUNCTION fkmin(use_pbc) bind (C,NAME='fkmin_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: use_pbc
 
if (use_pbc .eq. 0) then
  fkmin = 0.0
else
  fkmin = 0.0
  do l=1, ncells
    fkmin = fkmin + the_box(l)%minr
  enddo
  fkmin = fkmin / ncells
endif
 
END FUNCTION
 
REAL (KIND=c_double) FUNCTION oglmax() bind (C,NAME='oglmax_')
 
USE parameters
 
IMPLICIT NONE
 
DOUBLE PRECISION :: dmax
 
do l=1,3
  pmin(l) = fullpos(1,l,1)
  pmax(l) = pmin(l)
enddo
do i=1, ns
  if (i .eq. 1) then
    k = 2
  else
    k = 1
  endif
  do j=k, na
    do l=1,3
      pmin(l) = min(pmin(l), fullpos(j,l,i))
      pmax(l) = max(pmax(l), fullpos(j,l,i))
    enddo
  enddo
enddo
 
allocate (nfullpos(2,3,1))
do l=1, 3
  nfullpos(1,l,1) = pmin(l)
  nfullpos(2,l,1) = pmax(l)
enddo
call calcrij (1, 2, -2, 1, 1)
deallocate (nfullpos)
 
if (pbc) then
  oglmax = sqrt(dij)
  dmax = 0.0d0
  do l=1, ncells
    dmax = max(dmax, the_box(l)%maxv)
  enddo
  if (oglmax < dmax) oglmax = dmax*3.0;
else
  oglmax=2.0*sqrt(dij)
endif
 
if (oglmax .lt. 10.0) oglmax = 10.0
 
END FUNCTION