angles.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 bond_angles(nda) bind (C,NAME='bond_angles_')
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: nda
INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: anglea
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: sum_anga
DOUBLE PRECISION :: ang
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: angtab
#ifdef OPENMP
INTEGER :: numth
LOGICAL :: doatoms
#endif
INTERFACE
  DOUBLE PRECISION FUNCTION angijk (ATG1, ATG2, ATG3, ASTEP)
    INTEGER, INTENT(IN) :: atg1, atg2, atg3, astep
  END FUNCTION
END INTERFACE
 
if (allocated(anglea)) deallocate(anglea)
allocate(anglea(nsp,nsp,nsp,nda), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bond_angles"//char(0), "Table: ANGLEA"//char(0))
  bond_angles=0
  goto 001
endif
 
anglea(:,:,:,:)=0
 
delta_ang=180.0/dble(nda)
 
#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
  do i=1, ns
    ! OpemMP on Atoms
    !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
    !$OMP& PRIVATE(j, k, l, m, n, ANG, ANG_I) &
    !$OMP& SHARED(NUMTH, NS, i, NA, NCELLS, LOT, CONTJ, VOISJ, ANGLEA, DELTA_ANG, nda)
    !$OMP DO SCHEDULE(STATIC,NA/NUMTH)
    do j=1, na
 
      if (contj(j,i) .gt. 1) then
 
        do k=1, contj(j,i)-1
          do l=k+1, contj(j,i)
            m=voisj(k,j,i)
            n=voisj(l,j,i)
            ang = angijk(m, j, n, i)
 
            ang_i=anint(ang/delta_ang)
            if (ang_i.le.0) ang_i=1
            if (ang_i.gt.nda) ang_i=nda
            !$OMP ATOMIC
            anglea(lot(m),lot(j),lot(n),ang_i)=anglea(lot(m),lot(j),lot(n),ang_i)+1
            if (lot(m) .ne. lot(n)) then
              !$OMP ATOMIC
              anglea(lot(n),lot(j),lot(m),ang_i)=anglea(lot(n),lot(j),lot(m),ang_i)+1
            endif
          enddo
        enddo
 
      endif
 
    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, ANG, ANG_I) &
  !$OMP& SHARED(NUMTH, NS, NA, NCELLS, LOT, CONTJ, VOISJ, ANGLEA, DELTA_ANG, nda)
  !$OMP DO SCHEDULE(STATIC,NS/NUMTH)
#endif
  do i=1, ns
 
    do j=1, na
 
      if (contj(j,i) .gt. 1) then
 
        do k=1, contj(j,i)-1
          do l=k+1, contj(j,i)
            m=voisj(k,j,i)
            n=voisj(l,j,i)
            ang = angijk(m, j, n, i)
 
            ang_i=anint(ang/delta_ang)
            if (ang_i.le.0) ang_i=1
            if (ang_i.gt.nda) ang_i=nda
#ifdef OPENMP
            !$OMP ATOMIC
#endif
            anglea(lot(m),lot(j),lot(n),ang_i)=anglea(lot(m),lot(j),lot(n),ang_i)+1
            if (lot(m) .ne. lot(n)) then
#ifdef OPENMP
              !$OMP ATOMIC
#endif
              anglea(lot(n),lot(j),lot(m),ang_i)=anglea(lot(n),lot(j),lot(m),ang_i)+1
            endif
          enddo
        enddo
 
      endif
 
    enddo
  enddo
#ifdef OPENMP
  !$OMP END DO NOWAIT
  !$OMP END PARALLEL
endif
#endif
 
if (allocated(sum_anga)) deallocate(sum_anga)
allocate (sum_anga(nsp,nsp,nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bond_angles"//char(0), "Table: SUM_ANGA"//char(0))
  bond_angles=0
  goto 001
endif
 
do i=1, nsp
  do j=1, nsp
    do k=1, nsp
      sum_anga(i,j,k) = 0
      do l=1, nda
        sum_anga(i,j,k) = sum_anga(i,j,k) + anglea(i,j,k,l)
      enddo
    enddo
  enddo
enddo
 
if (allocated(angtab)) deallocate(angtab)
allocate(angtab(nda), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bond_angles"//char(0), "Table: ANGTAB"//char(0))
  bond_angles=0
  goto 001
endif
angtab(:)=0.0
 
m=0
do i=1, nsp
  do j=1, nsp
    do k=1, nsp
      if (sum_anga(i,j,k) .ne. 0) then
        do l=1, nda
          angtab(l)=100.0*dble(anglea(i,j,k,l))/dble(sum_anga(i,j,k))
        enddo
        call save_curve (nda, angtab, m, idan)
      else
        call save_curve (0, angtab, m, idan)
      endif
      m=m+1
    enddo
  enddo
enddo
 
bond_angles=1
 
001 continue
 
if (allocated(anglea)) deallocate(anglea)
if (allocated(sum_anga)) deallocate(sum_anga)
if (allocated(angtab)) deallocate(angtab)
 
END FUNCTION bond_angles
 
INTEGER (KIND=c_int) FUNCTION bond_diedrals(nda) bind (C,NAME='bond_diedrals_')
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: nda
INTEGER, DIMENSION(:,:,:,:,:), ALLOCATABLE :: angled
INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: sum_angd
DOUBLE PRECISION :: ang
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: angtab
#ifdef OPENMP
INTEGER :: numth
LOGICAL :: doatoms
#endif
INTERFACE
  DOUBLE PRECISION FUNCTION diedre (DG1, DG2, DG3, DG4, DSTEP)
    INTEGER, INTENT(IN) :: dg1, dg2, dg3, dg4, dstep
  END FUNCTION
END INTERFACE
 
if (allocated(angled)) deallocate(angled)
allocate(angled(nsp,nsp,nsp,nsp,nda), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bond_angles"//char(0), "Table: ANGLED"//char(0))
  bond_diedrals=0
  goto 001
endif
 
angled(:,:,:,:,:)=0
 
delta_ang=180.0/dble(nda)
 
#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, k, l, m, n, o, p, ANG, ANG_I) &
    !$OMP& SHARED(NUMTH, NS, i, NA, NCELLS, LOT, CONTJ, VOISJ, ANGLED, DELTA_ANG, nda)
    !$OMP DO SCHEDULE(STATIC,NA/NUMTH)
    do j=1, na
      do k=1, contj(j,i)
        m=voisj(k,j,i)
        if (contj(m,i) .ge. 2) then
          do l=1, contj(m,i)
            n=voisj(l,m,i)
            if (n .ne. j) then
              if (contj(n,i) .ge. 2) then
                do o=1, contj(n,i)
                  p = voisj(o,n,i)
                  if (p.ne.j .and. p.ne.m) then
                    ang=diedre(j, m, n, p, i)
                    ang_i=anint(ang/delta_ang)+1
                    if (ang_i.le.0) ang_i=1
                    if (ang_i.gt.nda) ang_i=nda
                    !$OMP ATOMIC
                    angled(lot(j),lot(m),lot(n),lot(p),ang_i)=angled(lot(j),lot(m),lot(n),lot(p),ang_i)+1
                    if (lot(j).ne.lot(m) .or. lot(j).ne.lot(n) .or. lot(j).ne.lot(p)) then
                     !$OMP ATOMIC
                      angled(lot(p),lot(n),lot(m),lot(j),ang_i)=angled(lot(p),lot(n),lot(m),lot(j),ang_i)+1
                    endif
                  endif
                enddo
              endif
            endif
          enddo
        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, ANG, ANG_I) &
  !$OMP& SHARED(NUMTH, NS, NA, NCELLS, LOT, CONTJ, VOISJ, ANGLED, DELTA_ANG, nda)
  !$OMP DO SCHEDULE(STATIC,NS/NUMTH)
#endif
  do i=1, ns
    do j=1, na
      do k=1, contj(j,i)
        m=voisj(k,j,i)
        if (contj(m,i) .ge. 2) then
          do l=1, contj(m,i)
            n=voisj(l,m,i)
            if (n .ne. j) then
              if (contj(n,i) .ge. 2) then
                do o=1, contj(n,i)
                  p = voisj(o,n,i)
                  if (p.ne.j .and. p.ne.m) then
                    ang=diedre(j, m, n, p, i)
                    ang_i=anint(ang/delta_ang)+1
                    if (ang_i.le.0) ang_i=1
                    if (ang_i.gt.nda) ang_i=nda
#ifdef OPENMP
                    !$OMP ATOMIC
#endif
                    angled(lot(j),lot(m),lot(n),lot(p),ang_i)=angled(lot(j),lot(m),lot(n),lot(p),ang_i)+1
                    if (lot(j).ne.lot(m) .or. lot(j).ne.lot(n) .or. lot(j).ne.lot(p)) then
#ifdef OPENMP
                     !$OMP ATOMIC
#endif
                      angled(lot(p),lot(n),lot(m),lot(j),ang_i)=angled(lot(p),lot(n),lot(m),lot(j),ang_i)+1
                    endif
                  endif
                enddo
              endif
            endif
          enddo
        endif
      enddo
    enddo
  enddo
#ifdef OPENMP
  !$OMP END DO NOWAIT
  !$OMP END PARALLEL
endif
#endif
if (allocated(sum_angd)) deallocate(sum_angd)
allocate (sum_angd(nsp,nsp,nsp,nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bond_angles"//char(0), "Table: SUM_ANGD"//char(0))
  bond_diedrals=0
  goto 001
endif
sum_angd(:,:,:,:) = 0
do i=1, nsp
  do j=1, nsp
    do k=1, nsp
      do l=1, nsp
        do m=1, nda
          sum_angd(i,j,k,l) = sum_angd(i,j,k,l) + angled(i,j,k,l,m)
        enddo
      enddo
    enddo
  enddo
enddo
 
if (allocated(angtab)) deallocate(angtab)
allocate(angtab(nda), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: bond_angles"//char(0), "Table: ANGTAB"//char(0))
  bond_diedrals=0
  goto 001
endif
angtab(:)=0.0
 
n=nsp*nsp*nsp
do i=1, nsp
  do j=1, nsp
    do k=1, nsp
      do l=1, nsp
        if (sum_angd(i,j,k,l) .ne. 0) then
          o=1
          do m=1, nda
            angtab(m)=100.0*dble(angled(i,j,k,l,m))/dble(sum_angd(i,j,k,l))
          enddo
          call save_curve (nda, angtab, n, idan)
        else
          call save_curve (0, angtab, n, idan)
        endif
        n=n+1
      enddo
    enddo
  enddo
enddo
 
bond_diedrals=1
 
001 continue
 
if (allocated(angled)) deallocate(angled)
if (allocated(sum_angd)) deallocate(sum_angd)
if (allocated(angtab)) deallocate(angtab)
 
END FUNCTION bond_diedrals