prepdata.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-2026 by CNRS and University of Strasbourg
!
!>
!! @file prepdata.F90
!! @short First level analysis of atomic coordinates
!! @author Sébastien Le Roux <sebastien.leroux@ipcms.unistra.fr>
 
#if defined (HAVE_CONFIG_H)
#  include <config.h>
#endif
 
SUBROUTINE free_contj_voisj () bind (C,NAME='free_contj_voisj_')
 
USE parameters
 
IMPLICIT NONE
 
if (allocated(voisj)) deallocate(voisj)
if (allocated(contj)) deallocate(contj)
 
END SUBROUTINE
 
SUBROUTINE read_contj (ATO, STP, CON) bind (C,NAME='read_contj_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: ATO, STP, CON
 
contj(ato+1,stp+1) = con
 
END SUBROUTINE
 
SUBROUTINE read_voisj (ATO, STP, CID, VID) bind (C,NAME='read_voisj_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: ATO, STP, CID, VID
 
voisj(cid+1,ato+1,stp+1) = vid+1
 
END SUBROUTINE
 
INTEGER (KIND=c_int) FUNCTION alloc_contj_voisj (N1, N2) bind (C,NAME='alloc_contj_voisj_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: n1, n2
 
! N1 = atomes, N2 = steps
if (allocated(voisj)) deallocate(voisj)
allocate(voisj(maxn,n1,n2), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: alloc_cont_vois"//char(0), "Table: VOISJ"//char(0))
  alloc_contj_voisj = 0
  goto 001
endif
if (allocated(contj)) deallocate(contj)
allocate(contj(n1,n2), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: alloc_cont_vois"//char(0), "Table: CONTJ"//char(0))
  alloc_contj_voisj = 0
endif
 
alloc_contj_voisj = 1
 
001 continue
 
END FUNCTION
 
INTEGER (KIND=c_int) FUNCTION alloc_data (N1, N2, N3) bind (C,NAME='alloc_data_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: n1, n2, n3
 
INTERFACE
  INTEGER FUNCTION allochem()
  END FUNCTION
END INTERFACE
 
na=n1
nsp=n2
ns=n3
 
if (allocated(fullpos)) deallocate(fullpos)
allocate (fullpos(na,3,ns), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: alloc_data"//char(0), "Table: FULLPOS"//char(0))
  alloc_data = 0
  goto 001
endif
 
if (allocated(lot)) deallocate(lot)
allocate(lot(na), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: alloc_data"//char(0), "Table: LOT"//char(0))
  alloc_data = 0
  goto 001
endif
alloc_data = allochem()
 
001 continue
 
END FUNCTION
 
SUBROUTINE prep_spec (idatoms, nsps, open_apf) bind (C,NAME='prep_spec_')
 
USE parameters
USE mendeleiev
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), DIMENSION(NSP), INTENT(IN) :: nsps
real(kind=c_double), DIMENSION(NSP), INTENT(IN) :: idatoms
INTEGER (KIND=c_int), INTENT(IN) :: open_apf
 
CHARACTER (LEN=14) :: ELEM
! Now we are calling the GTK+ routines
 
do i=1, nsp
  nbspbs(i) = nsps(i)
  j = int(idatoms(i))
  atomid(i) = j
  if (open_apf .eq. 1) then
    if (j .gt. 0) then
      tl(i) = atsym(j)
      elem = element(j)
    else
      tl(i) = "X "
      elem = "Unknown"
    endif
    ! In C all string must be terminated by a CHAR(0)
    ! To spec_data_
    call spec_data (0, i-1, atomid(i), nbspbs(i), &
                    tl(i)//char(0), elem//char(0), &
                    0.0d0, 0.0d0, 0.0d0, 0.0d0)
  endif
enddo
nbspbs(nsp+1) = na
 
END SUBROUTINE
 
SUBROUTINE read_chem (PMASS, PRAD, PNSCATT, PXSCATT) bind (C,NAME='read_chem_')
 
USE parameters
 
IMPLICIT NONE
 
real(kind=c_double), DIMENSION(NSP) :: pmass, prad, pnscatt, pxscatt
 
do i=1, nsp
  mass(i) = pmass(i)
  rvdw(i) = prad(i)
  nscattl(i) = pnscatt(i)
  xscattl(i) = pxscatt(i)
enddo
 
END SUBROUTINE
 
SUBROUTINE read_data (PLOT, PBS) bind (C,NAME='read_data_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), DIMENSION(NA), INTENT(IN) :: PLOT
INTEGER (KIND=c_int), DIMENSION(NSP), INTENT(IN) :: PBS
 
do i=1, na
  lot(i)=plot(i)+1
enddo
do i=1, nsp
  nbspbs(i) = pbs(i)
  xi(i) = dble(nbspbs(i)) / dble(na)
enddo
nbspbs(nsp+1) = na
 
END SUBROUTINE
 
SUBROUTINE read_pos (PCX, PCY, PCZ) bind (C,NAME='read_pos_')
 
USE parameters
 
real(kind=c_double), DIMENSION(NA*NS), INTENT(IN) :: pcx, pcy, pcz
 
k=0
do i=1, ns
  do j=1, na
    k=k+1
    fullpos(j,1,i)=pcx(k)
    fullpos(j,2,i)=pcy(k)
    fullpos(j,3,i)=pcz(k)
  enddo
enddo
 
END SUBROUTINE
 
INTEGER FUNCTION send_pos (NPA, NPS, NLOT, POSTAB)
 
INTEGER :: i, j, k, err
INTEGER, INTENT(IN) :: npa, nps
INTEGER, DIMENSION(:), INTENT(IN) :: nlot
DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN) :: postab
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xpos, ypos, zpos
 
if (allocated(xpos)) deallocate(xpos)
allocate(xpos(npa*nps), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: SEND_POS"//char(0), "Table: XPOS"//char(0))
  send_pos = 0
  goto 001
endif
if (allocated(ypos)) deallocate(ypos)
allocate(ypos(npa*nps), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: SEND_POS"//char(0), "Table: YPOS"//char(0))
  send_pos = 0
  goto 001
endif
if (allocated(zpos)) deallocate(zpos)
allocate(zpos(npa*nps), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: SEND_POS"//char(0), "Table: ZPOS"//char(0))
  send_pos = 0
  goto 001
endif
 
k=0
do i=1, nps
do j=1, npa
  k=k+1
  xpos(k)=postab(j,1,i)
  ypos(k)=postab(j,2,i)
  zpos(k)=postab(j,3,i)
enddo
enddo
 
! To 'save_pos_'
call save_pos (npa, nlot, k, xpos, ypos, zpos)
 
send_pos = 1
 
001 continue
 
if (allocated(xpos)) deallocate (xpos)
if (allocated(ypos)) deallocate (ypos)
if (allocated(zpos)) deallocate (zpos)
 
END FUNCTION
 
INTEGER (KIND=c_int) FUNCTION prep_data () bind (C,NAME='prep_data_')
 
!
! Data initialization
!
 
USE parameters
USE mendeleiev
 
IMPLICIT NONE
 
CHARACTER (LEN=2), DIMENSION(MAXE) :: ttype
 
INTERFACE
  INTEGER FUNCTION send_pos(NPA, NPS, NLOT, POSTAB)
    INTEGER, INTENT(IN) :: npa, nps
    INTEGER, DIMENSION(:), INTENT(IN) :: nlot
    DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN) :: postab
  END FUNCTION
  INTEGER FUNCTION allochem()
  END FUNCTION
  INTEGER FUNCTION findid(NAMEAT)
    CHARACTER (LEN=2), INTENT(IN) :: nameat
  END FUNCTION
  INTEGER FUNCTION chemistry ()
  END FUNCTION
END INTERFACE
 
prep_data = 0
 
if (allocated(lot)) deallocate(lot)
allocate(lot(na), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: prep_data"//char(0), "Table: LOT"//char(0))
  goto 001
endif
lot(:)=0        ! Initialisation du tableau
 
noc=1
do noa=1, na
  if (noa .eq. 1) then
    lot(noa)=noc
    ttype(noc)=tab_of_type(noa)
  else if (lot(noa) .eq. 0) then
    do nob=1, noc
      if (tab_of_type(noa) .eq. ttype(nob)) then
        lot(noa)= nob
      endif
    enddo
    if (lot(noa) .eq. 0) then
      noc=noc+1
      lot(noa)=noc
      ttype(noc)=tab_of_type(noa)
    endif
  endif
enddo
nsp=noc
 
if (allochem() == 0) then
  prep_data = 0
  goto 001
endif
 
if (allocated(elemid)) deallocate(elemid)
allocate(elemid(nsp), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: prep_data"//char(0), "Table: ELEMID"//char(0))
  goto 001
endif
 
nbspbs(:) = 0
do noa=1, na
  nob=lot(noa)
  nbspbs(nob)=nbspbs(nob)+1
enddo
 
nbspbs(nsp+1)=na
 
!
! Link between index and label
! Correspondance entre indice de type et label
!
 
do nob=1, nsp
  tl(nob)=ttype(nob)
enddo
 
do noa=1, nsp
  atomid(noa) = findid(tl(noa))
  select case (atomid(noa))
    case (-1)
      prep_data = 0
      goto 001
    case (0)
      elemid(noa) = "Dummy "//tl(noa)
      mass(noa) = 1.0
      rvdw(noa) = 0.5
      nscattl(noa) = 0.0
      xscattl(noa) = 0.0
    case default
      elemid(noa) = element(atomid(noa))
      mass(noa) = amass(atomid(noa))
      xscattl(noa) = atomid(noa)
      if(atomid(noa) < 105) then
        nscattl(noa) = coheb(atomid(noa))
        ! Covalent radius are the defaut values
        rvdw(noa) = arcov(atomid(noa))
        ! To use ionic radius uncomment the following line
        ! RVDW(NOA) = ARION(ATOMID(NOA))
        ! To use Wander Waals radius uncomment the following line
        ! RVDW(NOA) = ARVDW(ATOMID(NOA))
      else
        nscattl(noa) = 0.0d0
        rvdw(noa) = 0.0d0
      endif
  end select
  if (nscattl(noa) .eq. 0.00) then
      call show_warning ("Element "//tl(noa)//" does not have neutron scattering length "//char(0), &
                         "If this is a bug please report it to"//char(0), package_bugreport//char(0))
  endif
enddo
 
! Now we are calling the GTK+ routines
! To init_data_
call init_data (na, nsp, ns, 1)
 
do i=1, nsp
! In C all string must be terminated by a CHAR(0)
! To spec_data_
  call spec_data (1, i-1, atomid(i), nbspbs(i), &
                  tl(i)//char(0), elemid(i)//char(0), &
                  mass(i), rvdw(i), nscattl(i), xscattl(i))
enddo
 
if (chemistry() .eq. 1) then
  prep_data = send_pos(na, ns, lot, fullpos)
endif
 
001 continue
 
if (allocated(elemid)) deallocate(elemid)
 
END FUNCTION
 
INTEGER FUNCTION findid(NAMEAT)
 
USE parameters
USE mendeleiev
 
INTEGER (KIND=c_int), EXTERNAL :: dummy_ask
INTEGER :: elemt
CHARACTER (LEN=2), INTENT(IN) :: nameat
 
findid=0
do elemt=1, maxe
 
 if (atsym(elemt) .eq. nameat) then
   findid=elemt
   exit
 endif
 
enddo
 
if (findid .eq. 0) then
  call show_warning ("Problem with the atomic coordinates"//char(0), &
                     "Element "//nameat//" does not exist "//char(0), " "//char(0))
  findid = dummy_ask("Do you want to use dummy atom(s) for unknown species "//nameat//" ?"//char(0));
endif
 
END FUNCTION
 
REAL (KIND=c_double) FUNCTION set_mass (SW) bind (C,NAME='set_mass_')
 
USE mendeleiev
USE parameters
 
INTEGER (KIND=c_int), INTENT(IN) :: sw
 
set_mass = amass(sw)
 
END FUNCTION
 
REAL (KIND=c_double) FUNCTION set_radius (SW, RD) bind (C,NAME='set_radius_')
 
USE mendeleiev
USE parameters
 
INTEGER (KIND=c_int), INTENT(IN) :: sw, rd
DOUBLE PRECISION :: rad
 
if (sw .lt. 108) then
  if (rd .eq. 0) then
    rad = arcov(sw)
  else if (rd .eq. 1) then
    rad = arion(sw)
  else if (rd .eq. 2) then
    rad = arvdw(sw)
  else if (rd .eq. 3) then
    rad = arcry(sw)
  endif
else
  rad = 0.0
endif
set_radius = rad
 
END FUNCTION
 
REAL (KIND=c_double) FUNCTION set_neutron (SW) bind (C,NAME='set_neutron_')
 
USE mendeleiev
USE parameters
 
INTEGER (KIND=c_int), INTENT(IN) :: sw
 
if (sw .lt. 105) then
  set_neutron = coheb(sw)
else
  set_neutron = 0.0
endif
 
END FUNCTION
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                Reconstruction des trajectoires réelles
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE prep_pos (PINFO, FINFO) bind (C,NAME='prep_pos_')
 
USE parameters
 
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: PINFO, FINFO
 
INTERFACE
  INTEGER FUNCTION send_pos(NPA, NPS, NLOT, POSTAB)
    INTEGER, INTENT(IN) :: NPA, NPS
    INTEGER, DIMENSION(:), INTENT(IN) :: NLOT
    DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN) :: POSTAB
  END FUNCTION
END INTERFACE
 
pbc=.false.
if (pinfo .eq. 1) then
  pbc=.true.
endif
 
frac=.false.
if (finfo > 0) frac=.true.
 
if (frac) then
  do noc=1, ns
    do noa=1, na
      if (ncells .gt. 1) then
        fullpos(noa,:,noc) = matmul(fullpos(noa,:,noc),the_box(noc)%fractocart)
      else
        fullpos(noa,:,noc) = matmul(fullpos(noa,:,noc),the_box(1)%fractocart)
      endif
    enddo
  enddo
endif
 
noa = send_pos(na, ns, lot, fullpos)
 
END SUBROUTINE