molecules.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
!
 
LOGICAL FUNCTION add_mol(this_mol, MOL_ID, STEP_ID)
 
USE parameters
 
IMPLICIT NONE
 
TYPE (mol), POINTER, INTENT(INOUT) :: this_mol
INTEGER, INTENT(IN) :: mol_id, step_id
INTEGER :: aaa
 
TYPE (mol), POINTER :: new
 
allocate(new, stat=err)
if (err .ne. 0) then
  !call show_error ("Impossible to allocate memory"//CHAR(0), &
  !                 "Function: ADD_MOL"//CHAR(0), "Pointer: NEW"//CHAR(0))
  add_mol = .false.
  goto 001
endif
new%MID = mol_id
new%STEP = step_id
new%ATOMES = 0
allocate(new%BSP(nsp), stat=err)
if (err .ne. 0) then
  !call show_error ("Impossible to allocate memory"//CHAR(0), &
  !                 "Function: ADD_MOL"//CHAR(0), "TABLE: NEW%BSP"//CHAR(0))
  add_mol = .false.
  goto 001
endif
nullify(new%FIRST_AT)
do aaa=1, nsp
  new%BSP(aaa) = 0
enddo
nullify(new%NEXT)
new%PREV => this_mol
this_mol%NEXT => new
this_mol => new
add_mol = .true.
 
001 continue
 
END FUNCTION
 
LOGICAL FUNCTION add_ato(this_mol, ATOM_ID)
 
USE parameters
 
IMPLICIT NONE
 
TYPE (mol), POINTER, INTENT(INOUT) :: this_mol
INTEGER, INTENT(IN) :: atom_id
 
TYPE (at), POINTER :: new
 
allocate(new, stat=err)
if (err .ne. 0) then
  !call show_error ("Impossible to allocate memory"//CHAR(0), &
  !                 "Function: ADD_ATO"//CHAR(0), "POINTER: NEW"//CHAR(0))
  add_ato = .false.
  goto 001
endif
new%IND = atom_id
new%PREV => this_mol%ATOM
this_mol%ATOM%NEXT => new
this_mol%ATOM => new
add_ato=.true.
 
001 continue
 
END FUNCTION
 
RECURSIVE SUBROUTINE setmol (this_mol, toglin, stmb, molcount, the_step, the_atom, the_id)
 
USE parameters
 
IMPLICIT NONE
 
TYPE (mol), POINTER, INTENT(INOUT) :: this_mol
INTEGER, INTENT(INOUT) :: molcount, stmb
INTEGER, INTENT(IN) :: the_step, the_atom, the_id
INTEGER, DIMENSION(NA), INTENT(INOUT) :: toglin
INTEGER :: mc, md
INTERFACE
  LOGICAL FUNCTION add_ato(this_mol, ATOM_ID)
    USE parameters
    TYPE (mol), POINTER, INTENT(INOUT) :: this_mol
    INTEGER, INTENT(IN) :: atom_id
  END FUNCTION
END INTERFACE
 
molcount=molcount+1
toglin(the_atom) = the_id
this_mol%ATOMES = this_mol%ATOMES + 1
this_mol%BSP(lot(the_atom)) = this_mol%BSP(lot(the_atom)) + 1
stmb = stmb + 1
if (this_mol%ATOMES > 1) then
  if (.not. add_ato(this_mol, the_atom)) then
    molcount=-1
    goto 001
  endif
else
  allocate(this_mol%FIRST_AT, stat=err)
  if (err .ne. 0) then
    !call show_error ("Impossible to allocate memory"//CHAR(0), &
    !                 "Function: SETMOL"//CHAR(0), "POINTER: this_mol%FIRST_AT"//CHAR(0))
    molcount=-1
    goto 001
  endif
  !write (6, '("Adding first atom:: Mol= ",i4,", At= ",i4)') this_mol%MID, the_atom
  this_mol%FIRST_AT%IND = the_atom
  this_mol%ATOM => this_mol%FIRST_AT
endif
if  (stmb .eq. na) goto 002
if (molcount.eq.10000) goto 001
if (molcount.eq.-1) goto 001
do mc=1, contj(the_atom, the_step)
  md = voisj(mc,the_atom,the_step)
  if (toglin(md).eq.0 .and. molcount.lt.10000) then
    call setmol (this_mol, toglin, stmb, molcount, the_step, voisj(mc,the_atom,the_step), the_id)
  endif
  if (stmb .eq. na) goto 002
  if (molcount.eq.10000) goto 001
  if (molcount.eq.-1) goto 001
enddo
 
002 continue
 
molcount=molcount-1
 
001 continue
 
END SUBROUTINE setmol
 
INTEGER (KIND=c_int) FUNCTION molecules (frag_and_mol, allbonds) bind (C,NAME='molecules_')
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: frag_and_mol, allbonds
INTEGER :: totmol, molps, maxmol
INTEGER, DIMENSION(:), ALLOCATABLE :: atmol
INTEGER, DIMENSION(:), ALLOCATABLE :: atvs, mtmbs
#ifdef OPENMP
INTEGER :: numth
#endif
INTERFACE
  LOGICAL FUNCTION add_mol(this_mol, MOL_ID, STEP_ID)
    USE parameters
    TYPE (mol), POINTER, INTENT(INOUT) :: this_mol
    INTEGER, INTENT(IN) :: mol_id, step_id
  END FUNCTION
  RECURSIVE SUBROUTINE setmol (this_mol, toglin, stmb, molcount, the_step, the_atom, the_id)
    USE parameters
    TYPE (mol), POINTER, INTENT(INOUT) :: this_mol
    INTEGER, INTENT(INOUT) :: molcount, stmb
    INTEGER, INTENT(IN) :: the_step, the_atom, the_id
    INTEGER, DIMENSION(NA), INTENT(INOUT) :: toglin
  END SUBROUTINE
END INTERFACE
 
TYPE (at), POINTER :: tmpat
TYPE (mol), ALLOCATABLE, TARGET :: themol
TYPE (mol), POINTER :: tmpmol
 
if (allocated(fullpos)) deallocate(fullpos)
 
k = 0
do i=1, ns
 do j=1, na
   if (contj(j,i) .eq. 0) k = k + 1
 enddo
enddo
maxmol = k + allbonds/2;
 
if (allocated(themol)) deallocate(themol)
allocate(themol, stat=err)
if (err .ne. 0) then
  alc_tab="THEMOL"
  alc=.true.
  molecules = 0
  goto 001
endif
if (allocated(mtmbs)) deallocate(mtmbs)
allocate(mtmbs(ns), stat=err)
if (err .ne. 0) then
  alc_tab="MTMBS"
  alc=.true.
  molecules = 0
  goto 001
endif
 
if (frag_and_mol .eq. 1) then
  call allocate_mol_data ()
endif
 
molecules = 1
#ifdef OPENMP
numth = omp_get_max_threads()
if (ns.lt.numth) numth=ns
!$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
!$OMP& PRIVATE(i, j, k, l, m, n, o, ERR, TOGL, THEMOL, TMPMOL, &
!$OMP& TOTMOL, MOLCOUNTER, TMBS, ATMOL, TMPAT, ATVS) &
!$OMP& SHARED(NUMTH, frag_and_mol, NS, NA, NSP, LOT, MTMBS, CONTJ, VOISJ, ALC, ALC_TAB)
#endif
if (allocated(togl)) deallocate(togl)
allocate(togl(na), stat=err)
if (err .ne. 0) then
  alc_tab="TOGL"
  alc=.true.
  molecules = 0
#ifdef OPENMP
  goto 005
#else
  goto 001
#endif
endif
nullify(tmpat)
#ifdef OPENMP
  !$OMP DO SCHEDULE(STATIC,NS/NUMTH)
#endif
do i=1, ns
 
#ifdef OPENMP
  if (molecules .eq.0) goto 004
#endif
  togl(:)=0
  themol%MID = 1
  themol%STEP = i
  themol%ATOMES = 0
  nullify(themol%FIRST_AT)
  allocate(themol%BSP(nsp), stat=err)
  if (err .ne. 0) then
    alc_tab="THEMOL%BSP"
    alc=.true.
    molecules = 0
#ifdef OPENMP
    goto 004
#else
    goto 001
#endif
  endif
  do j=1, nsp
    themol%BSP(j) = 0
  enddo
  nullify(themol%NEXT)
  nullify(themol%PREV)
  tmbs=0
  totmol = 0
  mtmbs(i)=0
  molcounter=0
  tmpmol => themol
  n = 1
  002 continue
  do j=n, na
    if (molcounter.eq.0 .and. togl(j).eq.0) then
      totmol = totmol + 1
      if (totmol > 1) then
        if (.not.add_mol(tmpmol, mtmbs(i)+1, i)) then
          alc_tab="ADD_MOL"
          alc=.true.
          molecules = 0
#ifdef OPENMP
          goto 004
#else
          goto 001
#endif
        endif
      endif
      mtmbs(i) = mtmbs(i) + 1
      call setmol (tmpmol, togl, tmbs, molcounter, i, j, totmol)
      if (molcounter .lt. 0) then
        alc_tab="SETMOL"
        alc=.true.
#ifdef OPENMP
          goto 004
#else
          goto 001
#endif
      endif
      if (tmbs.eq.na) goto 003
      ! The 10000 iterations break required because of F90 limitations
      if (molcounter.eq.10000) goto 002
      molcounter = 0
      n = j
      goto 002
    else if (molcounter.gt.0 .and.togl(j).eq.totmol) then
      do l=1, contj(j,i)
        m = voisj(l,j,i)
        if (togl(m).eq.0) then
          molcounter = 0
          call setmol (tmpmol, togl, tmbs, molcounter, i, voisj(l,j,i), totmol)
          if (molcounter .lt. 0) then
            alc_tab="SETMOL"
            alc=.true.
#ifdef OPENMP
              goto 004
#else
              goto 001
#endif
          endif
          if (tmbs.eq.na) goto 003
          if (molcounter.eq.10000) goto 002
          molcounter = 1
        endif
      enddo
    endif
  enddo
 
  if (tmbs .lt. na) then
    molcounter = 0
    goto 002
  endif
 
  003 continue
 
  if (frag_and_mol .eq. 1) then
 
    tmpmol => themol
    call allocate_mol_for_step (i, mtmbs(i))
    do j=1, mtmbs(i)
      if (allocated(atmol)) deallocate(atmol)
      allocate(atmol(tmpmol%ATOMES), stat=err)
      if (err .ne. 0) then
        alc_tab="ATMOL"
        alc=.true.
        molecules = 0
#ifdef OPENMP
        goto 004
#else
        goto 001
#endif
      endif
      if (associated(tmpat)) deallocate (tmpat)
      allocate(tmpat, stat=err)
      if (err .ne. 0) then
        alc_tab="TMPAT"
        alc=.true.
        molecules = 0
#ifdef OPENMP
        goto 004
#else
        goto 001
#endif
      endif
      tmpat = tmpmol%FIRST_AT
      do k=1, tmpmol%ATOMES
        atmol(k) = tmpat%IND
        if (k < tmpmol%ATOMES) then
          tmpat => tmpat%NEXT
          deallocate (tmpat%PREV)
        else
          deallocate (tmpat)
        endif
      enddo
 
      call send_mol_details (tmpmol%STEP, tmpmol%MID, tmpmol%ATOMES, nsp, tmpmol%BSP, atmol)
      if (tmpmol%ATOMES .gt. 1) then
        do l=1, tmpmol%ATOMES
          m = atmol(l)
          n = contj(m,tmpmol%STEP)
          if (allocated(atvs)) deallocate(atvs)
          allocate(atvs(n), stat=err)
          if (err .ne. 0) then
            alc_tab="ATVS"
            alc=.true.
            molecules = 0
#ifdef OPENMP
            goto 004
#else
            goto 001
#endif
          endif
          do o=1, n
            atvs(o)= voisj(o,m,tmpmol%STEP)
          enddo
          call send_mol_neighbors (tmpmol%STEP, tmpmol%MID, m, n, atvs)
        enddo
      endif
      if (allocated(atvs)) deallocate(atvs)
      if (allocated(atmol)) deallocate(atmol)
      if (allocated(tmpmol%BSP)) deallocate (tmpmol%BSP)
 
      if (j .lt. mtmbs(i)) tmpmol => tmpmol%NEXT
 
    enddo
    call setup_molecules (i)
  endif
  call setup_fragments (i, togl)
  do while (tmpmol%MID .gt. 1)
    tmpmol => tmpmol%PREV
    deallocate (tmpmol%NEXT)
  enddo
#ifdef OPENMP
  004 continue
#endif
enddo
#ifdef OPENMP
!$OMP END DO NOWAIT
005 continue
if (allocated(togl)) deallocate (togl)
if (allocated(themol)) deallocate(themol)
if (allocated(atmol)) deallocate(atmol)
if (allocated(atvs)) deallocate(atvs)
!$OMP END PARALLEL
#else
if (allocated(togl)) deallocate (togl)
if (allocated(themol)) deallocate(themol)
if (allocated(atvs)) deallocate(atvs)
if (allocated(atmol)) deallocate(atmol)
#endif
 
molps = 0
j = 0
do i=1, ns
  molps = molps + mtmbs(i)
 j = max(mtmbs(i), j)
enddo
call send_coord_opengl (2, 1, 0, 0, j, 1)
call init_menu_fragmol (2)
 
if (frag_and_mol .eq. 1) then
  call setup_menu_molecules ()
endif
 
molecules=1
 
001 continue
 
if (alc) then
  call show_error ("Impossible to allocate memory !"//char(0), &
                   "Function: molecules"//char(0), char(9)//"Table: "//alc_tab(1:len_trim(alc_tab))//char(0))
endif
 
if (allocated(mtmbs)) deallocate(mtmbs)
 
END FUNCTION molecules