chemistry.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 chemistry () bind (C,NAME='chemistry_')
 
!
! Density, number density, concentrations, empirical formula
!
 
USE parameters
 
IMPLICIT NONE
 
INTEGER :: ref1, ref2, ref3
INTEGER, DIMENSION(:), ALLOCATABLE :: index_form           ! Formula index
INTEGER, DIMENSION(:), ALLOCATABLE :: refp
 
LOGICAL :: unclear, unfrac
 
!
! Determination of the empirical formula
! Détermination de la formule brute
!
if (allocated(index_form)) deallocate(index_form)
allocate(index_form(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: chemistry"//char(0), "Table: INDEX_FORM"//char(0))
  chemistry = 0
  goto 001
endif
 
if (allocated(refp)) deallocate(refp)
allocate(refp(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: chemistry"//char(0), "Table: REFP"//char(0))
  chemistry = 0
  goto 001
endif
 
do i=1, nsp
  refp(i)= nbspbs(i)
enddo
 
do i=1, nsp
 
  if (i .eq. 1)then
 
    ref1=nbspbs(i)
    ref2=i
 
  else
 
    ref3=ref1
    ref1=min(nbspbs(i),ref1)
    if (ref3 .ne. ref1) ref2=i
 
  endif
 
enddo
 
unclear=.true.
unfrac=.false.
l=1
do while (unclear)
 
  unclear=.false.
  do i=1, nsp
    if (i .ne. ref2) then
      if (.not.unfrac) then
        m=1
        do while (mod(m*refp(i),ref1).ne.0 .and. m.lt.ref1)
          m=m+1
        enddo
        if (mod(m*refp(i),ref1) .eq. 0) then
          index_form(i)=m*refp(i)/ref1
          unclear=.true.
        else
          index_form(i)=nbspbs(i)
          unclear=.true.
          unfrac=.true.
        endif
      else
        index_form(i)=nbspbs(i)
      endif
    endif
  enddo
  if (unfrac) then
    index_form(ref2)=nbspbs(ref2)
  else
    index_form(ref2)=refp(ref2)/nbspbs(ref2)
  endif
  if (unclear) then
    unclear=.false.
    do i=1, nsp-1
      do j=i+1, nsp
        y=dble(index_form(i))/dble(index_form(j))
        z=dble(nbspbs(i))/dble(nbspbs(j))
        if (y .ne. z) unclear=.true.
      enddo
    enddo
    if (unclear) then
      refp=refp/l
      l=l+1
      refp=refp*l
    endif
  endif
enddo
 
! To send_chem_info_
call send_chem_info (index_form)
 
chemistry = 1
 
001 continue
 
if (allocated(index_form)) deallocate(index_form)
if (allocated(refp)) deallocate(refp)
 
END FUNCTION