TABLE OF CONTENTS


ABINIT/m_ingeo [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ingeo

FUNCTION

 Initialize geometry variables for the ABINIT code.

COPYRIGHT

  Copyright (C) 1998-2022 ABINIT group (XG, RC)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_ingeo
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_atomdata
28  use m_sort
29  use m_dtset
30 
31  use m_symtk,      only : mati3inv, chkorthsy, symrelrot, mati3det, chkprimit, &
32 &                         symmetrize_rprimd, symmetrize_tnons,symmetrize_xred, symatm
33  use m_spgbuilder, only : gensymspgr, gensymshub, gensymshub4
34  use m_symfind,    only : symfind, symanal, symlatt
35  use m_geometry,   only : mkradim, mkrdim, xcart2xred, xred2xcart, randomcellpos, metric, reduce2primitive
36  use m_parser,     only : intagm, intagm_img, geo_t, geo_from_abivar_string, get_acell_rprim
37 
38  implicit none
39 
40  private

m_ingeo/fillcell [ Functions ]

[ Top ] [ m_ingeo ] [ Functions ]

NAME

 fillcell

FUNCTION

 Computes the atomic position of all the atoms in the unit cell starting
 with the symmetry operations and the atoms from the asymetric unit cell.

INPUTS

  chrgat(natom)=target charge for each atom. Not always used, it depends on the value of constraint_kind
  natrd = number of atoms in the assymetric unit cell
  natom = total number of atoms (to be checked)
  nsym = number of symmetry operations
  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
  symrel(3,3,nsym)=symmetry operations in real space in terms
   of primitive translations
  tnons(3,nsym)=nonsymmorphic translations for symmetry operations
  tolsym=tolerance on symmetries
  typat(1:natrd)=type integer for each atom in cell
  xred(3,1:natrd)=reduced dimensionless atomic coordinates

OUTPUT

SIDE EFFECTS

  At input, for the assymetric unit cell
  nucdipmom(3,1:natrd)=nuclear magnetic dipole moments of the atoms
  spinat(3,1:natrd)=spin-magnetization of the atoms
  typat(1:natrd)=type integer for each atom in cell
  xred(3,1:natrd)=reduced dimensionless atomic coordinates

  At output, for the complete unit cell
  nucdipmom(3,1:natom)=nuclear magnetic dipole moments of the atoms
  spinat(3,1:natom)=spin-magnetization of the atoms
  typat(1:natom)=type integer for each atom in cell
  xred(3,1:natom)=reduced dimensionless atomic coordinates

SOURCE

1909 subroutine fillcell(chrgat,natom,natrd,nsym,nucdipmom,spinat,symafm,symrel,tnons,tolsym,typat,xred)
1910 
1911 !Arguments ------------------------------------
1912 !scalars
1913  integer,intent(in) :: natom,natrd,nsym
1914 !arrays
1915  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
1916  integer,intent(inout) :: typat(natom)
1917  real(dp),intent(in) :: tolsym
1918  real(dp),intent(in) :: tnons(3,nsym)
1919  real(dp),intent(inout) :: chrgat(natom),nucdipmom(3,natom),spinat(3,natom),xred(3,natom)
1920 
1921 !Local variables ------------------------------
1922 !scalars
1923  integer :: curat,flagch,flageq,ii,iij,jj,kk
1924  character(len=500) :: msg
1925 !arrays
1926  integer :: bcktypat(nsym*natrd)
1927  real(dp) :: bckat(3),bcknucdipmom(3,nsym*natrd)
1928  real(dp) :: bckchrgat(nsym*natrd),bckspinat(3,nsym*natrd),bckxred(3,nsym*natrd)
1929 
1930 ! *************************************************************************
1931 
1932 !DEBUG
1933 !write(std_out,*)' fillcell : enter with nsym, natrd= ',nsym,natrd
1934 !write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
1935 !do ii=1,nsym
1936 !write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
1937 !end do
1938 !write(std_out,*)' Describe the input atoms (index,typat,xred,spinat)'
1939 !do jj=1,natrd
1940 !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj)
1941 !end do
1942 !ENDDEBUG
1943 
1944  curat=0
1945 
1946 !Cycle over all the symmetry operations
1947  do ii=1,nsym
1948 
1949 !  Cycle over all the atoms in the assymetric unit cell
1950    do jj=1,natrd
1951 
1952 !    Symmetry operation application
1953      bckat(:)=matmul(symrel(:,:,ii),xred(:,jj))+tnons(:,ii)
1954 
1955 !    Normalization of the coordinates in [0,1)
1956      do iij=1,3
1957        do while (bckat(iij)<-tolsym)
1958          bckat(iij)=bckat(iij)+1.0d0
1959        end do
1960        do while (bckat(iij)>=1.0d0-tolsym)
1961          bckat(iij)=bckat(iij)-1.0d0
1962        end do
1963      end do
1964 
1965 !    Check for duplicate atoms
1966      flagch=0
1967      do kk=1,curat
1968        flageq=0
1969        if ( abs(bckxred(1,kk)-bckat(1))<tolsym  .and. &
1970 &       abs(bckxred(2,kk)-bckat(2))<tolsym  .and. &
1971 &       abs(bckxred(3,kk)-bckat(3))<tolsym       ) exit
1972        flagch=flagch+1
1973      end do
1974 
1975      if (flagch==curat) then
1976 !      Add the obtained atom to the bckxred list
1977        curat=curat+1
1978        bckxred(:,curat)=bckat
1979        bcktypat(curat)=typat(jj)
1980        bckchrgat(curat)=chrgat(jj)
1981        bcknucdipmom(:,curat)=nucdipmom(:,jj)
1982        bckspinat(:,curat)=spinat(:,jj)*symafm(ii)
1983      end if
1984 
1985    end do
1986  end do
1987 
1988 !DEBUG
1989 !write(std_out,*)' fillcell : Proposed coordinates ='
1990 !do ii=1,curat
1991 !write(std_out,'(i4,3es16.6)' )ii,bckxred(:,ii)
1992 !end do
1993 !ENDDEBUG
1994 
1995  if (curat>natom) then
1996    write(msg, '(a,i3,a,a,i7,a,a,a,a)' )&
1997 &   'The number of atoms obtained from symmetries, ',curat,ch10,&
1998 &   'is greater than the input number of atoms, natom=',natom,ch10,&
1999 &   'This is not allowed.',ch10,&
2000 &   'Action: modify natom or the symmetry data in the input file.'
2001    ABI_ERROR(msg)
2002  end if
2003 
2004  if (curat<natom) then
2005    write(msg, '(a,i3,a,a,i7,a,a,a,a)' )&
2006 &   'The number of atoms obtained from symmetries, ',curat,ch10,&
2007 &   'is lower than the input number of atoms, natom=',natom,ch10,&
2008 &   'This is not allowed.',ch10,&
2009 &   'Action: modify natom or the symmetry data in the input file.'
2010    ABI_ERROR(msg)
2011  end if
2012 
2013 !Assignment of symmetry to xred
2014  xred(:,1:natom)=bckxred(:,1:natom)
2015  typat(1:natom)=bcktypat(1:natom)
2016  chrgat(1:natom)=bckchrgat(1:natom)
2017  nucdipmom(1:3,1:natom)=bcknucdipmom(1:3,1:natom)
2018  spinat(1:3,1:natom)=bckspinat(1:3,1:natom)
2019 
2020 !DEBUG
2021 !write(std_out,*)' fillcell : exit with natom=',natom
2022 !write(std_out,*)' Describe the output atoms (index,typat,xred,spinat)'
2023 !do jj=1,natom
2024 !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj)
2025 !end do
2026 !ENDDEBUG
2027 
2028 end subroutine fillcell

m_ingeo/ingeo [ Functions ]

[ Top ] [ m_ingeo ] [ Functions ]

NAME

 ingeo

FUNCTION

 Initialize geometry variables for the ABINIT code.
 1) set up unit cell: acell, rprim and rprimd ; deduce Bravais lattice
 2) (removed)
 3) Set up the number of atoms (natrd) in the primitive set, to be read.
 4) Read the type of each atom in the primitive set
 5) Read coordinates for each atom in the primitive set
 6) Eventually read the symmetries
 7) Checks whether the geometry builder must be used,
    and call it if needed. Call eventually the symmetry builder and analyser
    Make the adequate transfers if the geometry
    builder is not needed.
 8) Initialize the fixing of atoms, the initial velocities, and the initial atomic spin

INPUTS

 berryopt == 4/14: electric field is on; berryopt = 6/7/16/17: electric displacement field is on
 iimage= index of the current image
 iout=unit number of output file
 jdtset=number of the dataset looked for
 lenstr=actual length of the string
 msym=default maximal number of symmetries
 natom=number of atoms
 nimage=number of images
 npsp=number of pseudopotentials (needed for the dimension of znucl)
 nspden=number of spin-density components
 nsppol=number of independent spin polarizations
 ntypat=number of type of atoms
 nzchempot=defines the use of a spatially-varying chemical potential along z
 pawspnorb=1 when spin-orbit is activated within PAW
 ratsph(1:ntypat)=radius of the atomic sphere
 string*(*)=character string containing all the input data. Initialized previously in instrng.
 supercell_latt(3)=supercell lattice
 comm: MPI communicator

OUTPUT

 acell(3)=length of primitive vectors
 amu(ntypat)=mass of each atomic type
 bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
 chrgat(natom)=target charge for each atom. Not always used, it depends on the value of constraint_kind
 genafm(3)=magnetic translation generator (in case of Shubnikov group type IV)
 iatfix(3,natom)=indices for atoms fixed along some (or all) directions
 jellslab=not zero if jellslab keyword is activated
 slabzbeg, slabzend= the z coordinates of beginning / end of the jellium slab
 mixalch(npspalch,ntypalch)=alchemical mixing factors
 nsym=actual number of symmetries
 nucdipmom(3,natom)=nuclear magnetic dipole moment of each atom in atomic units
 ptgroupma = magnetic point group number
 rprim(3,3)=dimensionless real space primitive translations
 spgroup=symmetry space group
 spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
 symafm(1:msym)=(anti)ferromagnetic part of symmetry operations
 symmorphi=if 0, only allows symmorphic symmetry operations
 symrel(3,3,1:msym)=symmetry operations in real space in terms
  of primitive translations
 tnons(3,1:msym)=nonsymmorphic translations for symmetry operations
 tolsym=tolerance for the symmetry operations
 typat(natom)=type integer for each atom in cell
 vel(3,natom)=initial velocity of atoms in bohr/atomic time units
 vel_cell(3,3)=initial velocity of cell parameters in bohr/atomic time units
 xred(3,natom)=reduced dimensionless atomic coordinates
 znucl(1:npsp)=nuclear number of atom as specified in psp file

SIDE EFFECTS

NOTES

 the parameters ntypat and natom have already been read in indims,
 and were used to dimension the arrays needed here.

TODO

 The dtset datastructure should NOT be an argument of this routine ... !

 MG: I completely agree. Abinit developers must learn that Fortran does not allow for aliasing!

SOURCE

 130 subroutine ingeo (acell,amu,bravais,chrgat,dtset,&
 131   genafm,iatfix,icoulomb,iimage,iout,jdtset,jellslab,lenstr,mixalch,&
 132   msym,natom,nimage,npsp,npspalch,nspden,nsppol,nsym,ntypalch,ntypat,&
 133   nucdipmom,nzchempot,pawspnorb,&
 134   ptgroupma,ratsph,rprim,slabzbeg,slabzend,spgroup,spinat,string,supercell_lattice,symafm,&
 135   symmorphi,symrel,tnons,tolsym,typat,vel,vel_cell,xred,znucl,comm)
 136 
 137 !Arguments ------------------------------------
 138 !scalars
 139  integer,intent(in) :: iimage,iout,jdtset,lenstr,msym
 140  integer,intent(in) :: nimage,npsp,npspalch,nspden,nsppol
 141  integer,intent(in) :: ntypalch,ntypat,nzchempot,pawspnorb,comm
 142  integer,intent(inout) :: natom,symmorphi
 143  integer,intent(out) :: icoulomb,jellslab,ptgroupma,spgroup !vz_i
 144  integer,intent(inout) :: nsym !vz_i
 145  real(dp),intent(out) :: slabzbeg,slabzend,tolsym
 146  character(len=*),intent(in) :: string
 147 !arrays
 148  integer,intent(in) :: supercell_lattice(3)
 149  integer,intent(out) :: bravais(11),iatfix(3,natom) !vz_i
 150  integer,intent(inout) :: symafm(msym) !vz_i
 151  integer,intent(inout) :: symrel(3,3,msym) !vz_i
 152  integer,intent(out) :: typat(natom)
 153  real(dp),intent(inout) :: chrgat(natom)
 154  real(dp),intent(inout) :: nucdipmom(3,natom)
 155  real(dp),intent(in) :: ratsph(ntypat)
 156  real(dp),intent(inout) :: spinat(3,natom)
 157  real(dp),intent(out) :: acell(3),amu(ntypat),genafm(3),mixalch(npspalch,ntypalch)
 158  real(dp),intent(inout) :: rprim(3,3),tnons(3,msym) !vz_i
 159  real(dp),intent(out) :: vel(3,natom),vel_cell(3,3),xred(3,natom)
 160  real(dp),intent(in) :: znucl(npsp)
 161  type(dataset_type),intent(in) :: dtset
 162 
 163 !Local variables-------------------------------
 164  character(len=*), parameter :: format01110 ="(1x,a6,1x,(t9,8i8) )"
 165  character(len=*), parameter :: format01160 ="(1x,a6,1x,1p,(t9,3g18.10)) "
 166 !scalars
 167  integer, save :: print_comment_tolsym=1
 168  integer :: bckbrvltt,brvltt,chkprim,chkprim_fake,expert_user
 169  integer :: fixed_mismatch,i1,i2,i3,iatom,iatom_supercell,idir,ierr,iexit,ii
 170  integer :: ipsp,irreducible,isym,itranslat,itypat,jsym,marr,mismatch_fft_tnons,multi,multiplicity,natom_uc,natfix,natrd
 171  integer :: nobj,noncoll,nptsym,nsym_now,ntranslat,ntyppure,random_atpos,shubnikov,spgaxor,spgorig
 172  integer :: spgroupma,tgenafm,tnatrd,tread,try_primitive,tscalecart,tspgroupma, tread_geo
 173  integer :: txcart,txred,txrandom,use_inversion
 174  real(dp) :: amu_default,ucvol,sumalch
 175  character(len=1000) :: msg
 176  character(len=lenstr) :: geo_string
 177  type(atomdata_t) :: atom
 178  type(geo_t) :: geo
 179 !arrays
 180  integer :: bravais_reduced(11)
 181  integer,allocatable :: indsym(:,:,:),intarr(:)
 182  integer,allocatable :: is_translation(:)
 183  integer,allocatable :: ptsymrel(:,:,:),typat_read(:),symrec(:,:,:)
 184  real(dp) :: angdeg(3), field_xred(3),gmet(3,3),gprimd(3,3),rmet(3,3),rcm(3)
 185  real(dp) :: rprimd(3,3),rprimd_read(3,3),rprimd_new(3,3),rprimd_primitive(3,3),scalecart(3)
 186  real(dp),allocatable :: mass_psp(:)
 187  real(dp),allocatable :: tnons_cart(:,:),tnons_new(:,:),translations(:,:)
 188  real(dp),allocatable :: xcart(:,:),xcart_read(:,:),xred_read(:,:),dprarr(:)
 189 
 190 ! *************************************************************************
 191 
 192 !DEBUG
 193 !write(std_out,'(a)')' m_ingeo%ingeo : enter '
 194 !call flush(std_out)
 195 !ENDDEBUG
 196 
 197  marr=max(12,3*natom,9*msym)
 198  ABI_MALLOC(intarr,(marr))
 199  ABI_MALLOC(dprarr,(marr))
 200 
 201  ! Try from geo_string
 202  call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'structure', tread_geo, 'KEY', key_value=geo_string)
 203 
 204  if (tread_geo /= 0) then
 205    ! Set up unit cell from external file.
 206    geo = geo_from_abivar_string(geo_string, comm)
 207    acell = one
 208    rprim = geo%rprimd
 209    !call exclude(lenstr, string, jdtset, iimage, "acell, rprim, angdeg, scalecar")
 210 
 211  else
 212    ! Set up unit cell from acell, rprim, angdeg
 213    call get_acell_rprim(lenstr, string, jdtset, iimage, nimage, marr, acell, rprim)
 214  end if ! geo% or (acell, rprim, angdeg)
 215 
 216  ! Rescale rprim using scalecart (and set scalecart to one)
 217  scalecart(1:3)=one
 218  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'scalecart',tscalecart,'LEN')
 219  if(tscalecart==1) scalecart(1:3)=dprarr(1:3)
 220  call intagm_img(scalecart,iimage,jdtset,lenstr,nimage,3,string,"scalecart",tscalecart,'LEN')
 221 
 222  rprim(:,1)=scalecart(:)*rprim(:,1)
 223  rprim(:,2)=scalecart(:)*rprim(:,2)
 224  rprim(:,3)=scalecart(:)*rprim(:,3)
 225  scalecart(:)=one
 226 
 227  ! Compute the multiplicity of the supercell
 228  multiplicity=supercell_lattice(1)*supercell_lattice(2)*supercell_lattice(3)
 229 
 230  if (tread_geo == 0) then
 231    ! Get the number of atom in the unit cell. Read natom from string
 232    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natom',tread,'INT')
 233 
 234    ! Might initialize natom from XYZ file
 235    if (tread==0) call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'_natom',tread,'INT')
 236 
 237    if(tread==1) natom_uc=intarr(1)
 238  else
 239    natom_uc = geo%natom
 240  end if
 241 
 242  ! Store the rprimd of the unit cell
 243  call mkrdim(acell,rprim,rprimd_read)
 244 
 245  ! Multiply the rprim to get the rprim of the supercell
 246  if(multiplicity > 1)then
 247    rprim(:,1) = rprim(:,1) * supercell_lattice(1)
 248    rprim(:,2) = rprim(:,2) * supercell_lattice(2)
 249    rprim(:,3) = rprim(:,3) * supercell_lattice(3)
 250  end if
 251 
 252  ! Compute different matrices in real and reciprocal space, also checks whether ucvol is positive.
 253  call mkrdim(acell, rprim, rprimd)
 254  call metric(gmet, gprimd, -1, rmet, rprimd, ucvol)
 255 
 256 !tolsym = tol8
 257 !XG20200801 New default value for tolsym. This default value is also defined in m_invars1.F90
 258  tolsym = tol5
 259  !if (tread_geo /= 0 .and. geo%filetype == "poscar") then
 260  !  tolsym = tol4
 261  !  ABI_COMMENT("Reading structure from POSCAR --> default value of tolsym is set to 1e-4")
 262  !end if
 263 
 264  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'tolsym',tread,'DPR')
 265  if(tread==1) tolsym=dprarr(1)
 266 
 267  ! Find a tentative Bravais lattice and its point symmetries (might not use them)
 268  ! Note that the Bravais lattice might not be the correct one yet (because the
 269  ! actual atomic locations might lower the symmetry obtained from the lattice parameters only)
 270  ABI_MALLOC(ptsymrel,(3,3,msym))
 271  call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
 272 
 273  ! 3) Possibly, initialize a jellium slab
 274  jellslab=0
 275  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'jellslab',tread,'INT')
 276  if(tread==1) jellslab=intarr(1)
 277 
 278  slabzbeg=zero
 279  slabzend=zero
 280  if(jellslab/=0)then
 281    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'slabzbeg',tread,'DPR')
 282    if(tread==1) slabzbeg=dprarr(1)
 283 
 284    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'slabzend',tread,'DPR')
 285    if(tread==1) slabzend=dprarr(1)
 286  end if
 287 
 288  ! 4) Set up the number of atoms in the primitive set, to be read.
 289  ! This is the default
 290  natrd=natom
 291  if(multiplicity > 1) natrd = natom_uc
 292 
 293  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natrd',tnatrd,'INT')
 294  if(tnatrd==1) natrd=intarr(1)
 295 
 296  if(natrd<1 .or. natrd>natom)then
 297    if(natrd>1 .and. multiplicity > 1) then
 298      if(natrd < natom)then
 299        write(msg, '(3a)' )&
 300         'The number of atoms to be read (natrd) can not be used with supercell_latt.',ch10,&
 301         'Action: Remove natrd or supercell_latt in your input file.'
 302        ABI_ERROR(msg)
 303      else
 304        write(msg,'(3a,I0,a,I0,a,I0,2a)')&
 305        'The input variable supercell_latt is present',ch10,&
 306        'thus a supercell ',supercell_lattice(1),' ',supercell_lattice(2),&
 307        ' ',supercell_lattice(3),' is generated',ch10
 308        ABI_WARNING(msg)
 309      end if
 310    else
 311      write(msg, '(3a,i0,a,i0,2a,a)' )&
 312       'The number of atoms to be read (natrd) must be positive and not bigger than natom.',ch10,&
 313       'This is not the case: natrd=',natrd,', natom=',natom,ch10,&
 314       'Action: correct natrd or natom in your input file.'
 315      ABI_ERROR(msg)
 316    end if
 317  end if
 318 
 319  ! 5) Read the type and initial spin of each atom in the primitive set--------
 320  ABI_MALLOC(typat_read,(natrd))
 321  typat_read(1)=1
 322 
 323  if (tread_geo == 0) then
 324    call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'typat',tread,'INT')
 325 
 326    ! If not read, try the XYZ data
 327    if(tread==0) call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'_typat',tread,'INT')
 328    if(tread==1) typat_read(1:natrd)=intarr(1:natrd)
 329 
 330  else
 331    typat_read = geo%typat
 332  end if
 333 
 334  do iatom=1,natrd
 335    if(typat_read(iatom)<1 .or. typat_read(iatom)>ntypat )then
 336      write(msg,'(a,i0,a,i0,a,a,a,i0,a,a,a)')&
 337       'The input type of atom number ',iatom,' is equal to ',typat_read(iatom),',',ch10,&
 338       'while it should be between 1 and ntypat= ',ntypat,'.',ch10,&
 339       'Action: change either the variable typat or the variable ntypat.'
 340      ABI_ERROR(msg)
 341    end if
 342  end do
 343 
 344  ! 6) Read coordinates for each atom in the primitive set--------
 345 
 346  ABI_MALLOC(xcart_read,(3,natrd))
 347  ABI_MALLOC(xred_read,(3,natrd))
 348 
 349  random_atpos=0
 350  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'random_atpos',txrandom,'INT')
 351  if(txrandom==1) random_atpos=intarr(1)
 352  if (random_atpos < 0 .or. random_atpos > 5) then
 353    write(msg,'(3a)')&
 354     'Random positions is a variable defined between 0 and 5. Error in the input file. ',ch10,&
 355     'Action: define one of these in your input file.'
 356    ABI_ERROR(msg)
 357  end if
 358 
 359  !if(nimage/=1 .and. iimage/=1)then
 360  !FIXME: should this be called outside the above end if?
 361  call randomcellpos(natom,npsp,ntypat,random_atpos,ratsph,rprim,rprimd_read,typat_read,&
 362                     xred_read(:,1:natrd),znucl,acell)
 363  !This should not be printed if randomcellpos did nothing - it contains garbage. Spurious output anyway
 364  !end if
 365 
 366  if (tread_geo == 0) then
 367 
 368    call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xred',txred,'DPR')
 369    if (txred==1 .and. txrandom == 0) xred_read(:,1:natrd) = reshape(dprarr(1:3*natrd) , [3, natrd])
 370    call intagm_img(xred_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xred",txred,'DPR')
 371 
 372    call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xcart',txcart,'LEN')
 373    if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd])
 374    call intagm_img(xcart_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xcart",txcart,'LEN')
 375 
 376    ! Might initialize xred from XYZ file
 377    if (txred+txcart+txrandom==0) then
 378      call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'_xred',txred,'DPR')
 379      if (txred==1 .and. txrandom==0) xred_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd])
 380 
 381      call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'_xcart',txcart,'DPR')
 382      if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd])
 383    end if
 384 
 385  else
 386    txcart = 0; txrandom = 0; txred = 1
 387    xred_read = geo%xred
 388  end if
 389 
 390  if (txred + txcart + txrandom == 0) then
 391    write(msg, '(3a)' )&
 392     'Neither xred nor xcart are present in input file. ',ch10,&
 393     'Action: define one of these in your input file.'
 394    ABI_ERROR(msg)
 395  end if
 396 
 397  if (txred==1)   write(msg, '(a)' ) '  xred   is defined in input file'
 398  if (txcart ==1) write(msg, '(a)' ) '  xcart  is defined in input file (possibly in Angstrom)'
 399  if (txrandom ==1) write(msg, '(a)' ) '  xred  as random positions in the unit cell'
 400  if (txrandom ==1) write(msg, '(a)' ) '  xcart  are defined from a random distribution '
 401  call wrtout(std_out, msg)
 402 
 403  if (txred + txcart + txrandom > 1)then
 404    write(msg, '(3a)' )&
 405     'Too many input channels for atomic positions are defined.',ch10,&
 406     'Action: choose to define only one of these.'
 407    ABI_ERROR(msg)
 408  end if
 409 
 410  if (txred==1 .or. txrandom /=0 ) then
 411    call wrtout(std_out,' ingeo: takes atomic coordinates from input array xred ')
 412    call xred2xcart(natrd,rprimd_read,xcart_read,xred_read)
 413  else
 414    call wrtout(std_out,' ingeo: takes atomic coordinates from input array xcart')
 415    txcart=1
 416  end if
 417 
 418  !At this stage, the cartesian coordinates are known, for the atoms whose coordinates where read.
 419 
 420  ! Here, allocate the variable that will contain the completed
 421  ! sets of xcart, after the use of the geometry builder or the symmetry builder
 422  ABI_MALLOC(xcart,(3,natom))
 423 
 424  !7) Eventually read the symmetries
 425  !Take care of the symmetries
 426 
 427  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nsym',tread,'INT')
 428  if(tread==1) nsym=intarr(1)
 429 
 430  ! Check that nsym is not negative
 431  if (nsym<0) then
 432    write(msg, '(a,i0,4a)' )&
 433     'Input nsym must be positive or 0, but was ',nsym,ch10,&
 434     'This is not allowed.',ch10,'Action: correct nsym in your input file.'
 435    ABI_ERROR(msg)
 436  end if
 437  ! Check that nsym is not bigger than msym
 438  if (nsym>msym) then
 439    write(msg, '(2(a,i0),5a)')&
 440     'Input nsym = ',nsym,' exceeds msym = ',msym,'.',ch10,&
 441     'This is not allowed.',ch10,'Action: correct nsym in your input file.'
 442    ABI_ERROR(msg)
 443  end if
 444  if (multiplicity>1) then
 445    nsym = 1
 446    ABI_WARNING('Input nsym is now set to one due to the supercell_latt input')
 447  end if
 448 
 449  ! Read symmorphi
 450  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'symmorphi',tread,'INT')
 451  if(tread==1) symmorphi=intarr(1)
 452 
 453  ! Now, read the symmetry operations
 454  if(nsym>0)then
 455    call intagm(dprarr,intarr,jdtset,marr,9*nsym,string(1:lenstr),'symrel',tread,'INT')
 456    if(nsym>1 .and. tread==0)then
 457      write(msg,'(3a)')&
 458        'When nsym>1, symrel must be defined in the input file.',ch10,&
 459        'Action: either change nsym, or define symrel in your input file.'
 460      ABI_ERROR(msg)
 461    end if
 462    if(tread==1) symrel(:,:,1:nsym)=reshape( intarr(1:9*nsym) , [3, 3, nsym])
 463 
 464    ! Take care of tnons
 465    tnons(:,1:nsym)=zero
 466    call intagm(dprarr,intarr,jdtset,marr,3*nsym,string(1:lenstr),'tnons',tread,'DPR')
 467    if(tread==1) tnons(:,1:nsym)=reshape( dprarr(1:3*nsym), [3, nsym])
 468 
 469    if(symmorphi==0)then
 470      do isym=1,nsym
 471        if(sum(tnons(:,isym)**2)>tol6)then
 472          write(msg, '(5a,i0,a,3f8.4,3a)' )&
 473          'When symmorphi /= 1, the vectors of translation (tnons)',ch10,&
 474          'a symmetry operation must vanish.',ch10,&
 475          'However, for the symmetry operation number ',isym,', tnons =',tnons(:,isym),'.',ch10,&
 476          'Action: either change your list of allowed symmetry operations, or use the symmetry finder (nsym=0).'
 477          ABI_ERROR(msg)
 478        end if
 479      end do
 480    end if
 481 
 482    ! Take care of symafm
 483    call intagm(dprarr,intarr,jdtset,marr,nsym,string(1:lenstr),'symafm',tread,'INT')
 484    if(tread==1) symafm(1:nsym)=intarr(1:nsym)
 485  end if
 486 
 487 
 488  !8) Checks whether the geometry builder must be used, and call it if needed.
 489  !Call the symmetry builder and analyzer if needed.
 490 
 491  ! At this stage, nsym might still contain the default 0, msym contains the default dtset%maxnsym.
 492  ! The cartesian coordinates of the atoms of the primitive set are contained in xcart_read.
 493 
 494  nobj=0
 495  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nobj',tread,'INT')
 496  if(tread==1) nobj=intarr(1)
 497  if(nobj /= 0 .and. multiplicity > 1)then
 498    write(msg, '(3a)' )&
 499     'nobj can not be used with supercell_latt.',ch10,&
 500     'Action: Remove nobj or supercell_latt in your input file.'
 501    ABI_ERROR(msg)
 502  end if
 503 
 504 !If there are objects, chkprim will not be used immediately
 505 !But, if there are no objects, but a space group, it will be used directly.
 506 !Need first to check the value of expert_user
 507  expert_user=0
 508  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'expert_user',tread,'INT')
 509  if(tread==1) expert_user=intarr(1)
 510  if(expert_user==0)then
 511    chkprim=1
 512    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'chkprim',tread,'INT')
 513    if(tread==1) chkprim=intarr(1)
 514  else
 515    chkprim=0
 516  endif
 517 
 518  if(nobj/=0)then
 519 
 520    ! chrgat is read for each atom, from 1 to natom
 521    call intagm(dprarr,intarr,jdtset,marr,natom,string(1:lenstr),'chrgat',tread,'DPR')
 522    if(tread==1) then
 523      chrgat(1:natom) = dprarr(1:natom)
 524    end if
 525 
 526    ! Spinat is read for each atom, from 1 to natom
 527    call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'spinat',tread,'DPR')
 528    if(tread==1) then
 529      spinat(1:3,1:natom) = reshape( dprarr(1:3*natom) , [3, natom])
 530    else if (nspden==4.or.(nspden==2.and.nsppol==1)) then
 531      write(msg, '(5a)' )&
 532       'When nspden=4 or (nspden==2 and nsppol==1), the input variable spinat must be',ch10,&
 533       'defined in the input file, which is apparently not the case.',ch10,&
 534       'Action: define spinat or use nspden=1 in your input file.'
 535      ABI_ERROR(msg)
 536    end if
 537 
 538    ! nucdipmom is read for each atom, from 1 to natom
 539    call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'nucdipmom',tread,'DPR')
 540    if(tread==1) then
 541      nucdipmom(1:3,1:natom) = reshape( dprarr(1:3*natom), [3, natom])
 542    end if
 543 
 544    ! Will use the geometry builder
 545    if(tnatrd/=1 .and. nobj/=0)then
 546      write(msg, '(3a,i0,5a)' )&
 547       'The number of atoms to be read (natrd) must be initialized',ch10,&
 548       'in the input file, when nobj= ',nobj,'.',ch10,&
 549       'This is not the case.',ch10,&
 550       'Action: initialize natrd in your input file.'
 551      ABI_ERROR(msg)
 552    end if
 553 
 554    if(jellslab/=0)then
 555      write(msg, '(a,i0,3a)' )&
 556       'A jellium slab cannot be used when nobj= ',nobj,'.',ch10,&
 557       'Action: change one of the input variables jellslab or nobj in your input file.'
 558      ABI_ERROR(msg)
 559    end if
 560 
 561    call ingeobld (iout,jdtset,lenstr,natrd,natom,nobj,string,typat,typat_read,xcart,xcart_read)
 562 
 563    ! Finalize the computation of coordinates: produce xred.
 564    call xcart2xred(natom,rprimd,xcart,xred)
 565 
 566  else
 567    ! nobj==0
 568 
 569    ! chrgat is read for each irreducible atom, from 1 to natrd
 570    call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'chrgat',tread,'DPR')
 571    if(tread==1)chrgat(1:natrd) = dprarr(1:natrd)
 572 
 573    ! Spinat is read for each irreducible atom, from 1 to natrd
 574    call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'spinat',tread,'DPR')
 575    if(tread==1)spinat(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natrd])
 576 
 577    ! nucdipmom is read for each irreducible atom, from 1 to natrd
 578    call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'nucdipmom',tread,'DPR')
 579    if(tread==1)nucdipmom(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natrd])
 580 
 581    ! Compute xred/typat and spinat for the supercell
 582    if(multiplicity > 1)then
 583      iatom_supercell = 0
 584      do i1 = 1, supercell_lattice(1)
 585        do i2 = 1, supercell_lattice(2)
 586          do i3 = 1, supercell_lattice(3)
 587            do iatom = 1, natom_uc
 588              iatom_supercell = iatom_supercell + 1
 589              xcart(:,iatom_supercell) = xcart_read(:,iatom) + matmul(rprimd_read,(/i1-1,i2-1,i3-1/))
 590              chrgat(iatom_supercell) = chrgat(iatom)
 591              spinat(1:3,iatom_supercell) = spinat(1:3,iatom)
 592              typat(iatom_supercell) = typat_read(iatom)
 593            end do
 594          end do
 595        end do
 596      end do
 597      call xcart2xred(natom,rprimd,xcart,xred)
 598    else
 599      ! No supercell
 600      call xcart2xred(natrd,rprimd,xcart_read,xred)
 601    end if
 602 
 603    spgroup=0
 604    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgroup',tread,'INT')
 605    if(tread==1) spgroup=intarr(1)
 606 
 607    if(spgroup/=0 .or. nsym/=0)then
 608 
 609      if(jellslab/=0 .and. nsym/=1 .and. spgroup/=1)then
 610        write(msg, '(5a)' )&
 611         'For the time being, a jellium slab can only be used',ch10,&
 612         'either with the symmetry finder (nsym=0) or with the space group 1 (nsym=1)',ch10,&
 613         'Action: change one of the input variables jellslab or nsym or spgroup in your input file.'
 614        ABI_ERROR(msg)
 615      end if
 616 
 617      if(nzchempot/=0 .and. nsym/=1 .and. spgroup/=1)then
 618        write(msg, '(5a)' )&
 619         'For the time being, a spatially-varying chemical potential can only be used',ch10,&
 620         'either with the symmetry finder (nsym=0) or with the space group 1 (nsym=1)',ch10,&
 621         'Action: change one of the input variables nzchempot or nsym or spgroup in your input file.'
 622        ABI_ERROR(msg)
 623      end if
 624 
 625      typat(1:natrd)=typat_read(1:natrd)
 626 
 627      if(spgroup/=0 .and. nsym/=0)then
 628        write(msg, '(a,i0,a,a,i0,a,a,a,a,a,a,a,a)' )&
 629          'The spatial group number spgroup= ',spgroup,ch10,&
 630          'is specified, as well as the number of symmetries nsym= ',nsym,ch10,&
 631          'This is not allowed, as you can define the symmetries',ch10,&
 632          'either using spgroup OR using nsym, but not both.',ch10,&
 633          'Action: modify your input file',ch10,&
 634          '(either set spgroup to 0, or nsym to 0)'
 635        ABI_ERROR(msg)
 636      end if
 637 
 638      brvltt=0
 639 
 640      if(spgroup/=0)then
 641 
 642        ! Will generate the spatial group using spgroup
 643        ! Assign default values
 644        spgaxor=1
 645        spgorig=1
 646        call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'brvltt',tread,'INT')
 647        if(tread==1) brvltt=intarr(1)
 648        call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgaxor',tread,'INT')
 649        if(tread==1) spgaxor=intarr(1)
 650        call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgorig',tread,'INT')
 651        if(tread==1) spgorig=intarr(1)
 652 
 653        ! Treat the case of magnetic groups
 654        call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgroupma',tspgroupma,'INT')
 655        if(tspgroupma==1) spgroupma=intarr(1)
 656        call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'genafm',tgenafm,'DPR')
 657        if(tgenafm==1) genafm(1:3)=dprarr(1:3)
 658        if(tspgroupma/=0 .and. tgenafm/=0)then
 659          write(msg, '(a,i0,a,a,3es9.2,a,a,a,a,a,a,a,a)' )&
 660            'The spatial group number spgroupma= ',spgroupma,ch10,&
 661            'is specified, as well as the antiferromagnetic generator genafm=',genafm(1:3),ch10,&
 662            'This is not allowed, as you can define the magnetic space group',ch10,&
 663            'either using spgroupma OR using genafm, but not both.',ch10,&
 664            'Action: modify your input file',ch10,&
 665            '(either define spgroupma or genafm)'
 666          ABI_ERROR(msg)
 667        end if
 668 
 669        ! TODO: all the symmetry generation operations should be in one big routine
 670 
 671        ! If spgroupma is defined, check whether it is consistent
 672        ! with spgroup, determine the Shubnikov type,
 673        ! and, for type IV, find the corresponding genafm
 674        shubnikov=1
 675        if(tspgroupma==1)then
 676          call gensymshub(genafm,spgroup,spgroupma,shubnikov)
 677        else if(tgenafm==1)then
 678          shubnikov=4
 679        end if
 680 
 681        ! Generate the spatial group of symmetries in a conventional cell
 682        ! In case of Shubnikov space group type IV, only generate the
 683        ! Fedorov (non-magnetic) group. For Shubnikov type III space group,
 684        ! the magnetic part is generated here.
 685        bckbrvltt=brvltt
 686        if(brvltt==-1)brvltt=0
 687        call gensymspgr(brvltt,msym,nsym,shubnikov,spgaxor,spgorig,spgroup,spgroupma,symafm,symrel,tnons)
 688 
 689        ! For shubnikov type IV groups,
 690        ! double the space group, using the antiferromagnetic translation generator
 691        if(shubnikov==4)then
 692          call gensymshub4(genafm,msym,nsym,symafm,symrel,tnons)
 693        end if
 694 
 695        !write(std_out,*)' after gensymshub4, nsym =',nsym
 696        !write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
 697        !do ii=1,nsym
 698        !write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
 699        !end do
 700 
 701        ! If brvltt was -1 at input, one should now change the conventional cell
 702        ! to a primitive one, if brvltt/=1
 703        if(bckbrvltt==-1 .and. brvltt/=1)then
 704          ! Will work with rprim only
 705          rprim(:,:)=rprimd(:,:)
 706          rprimd_new(:,:)=rprimd(:,:)
 707          acell(:)=1.0_dp
 708          select case(brvltt)
 709          case(5)
 710            rprimd_new(:,2)=(rprim(:,2)+rprim(:,3))*0.5_dp
 711            rprimd_new(:,3)=(rprim(:,3)-rprim(:,2))*0.5_dp
 712          case(6)
 713            rprimd_new(:,1)=(rprim(:,1)+rprim(:,3))*0.5_dp
 714            rprimd_new(:,3)=(rprim(:,3)-rprim(:,1))*0.5_dp
 715          case(4)
 716            rprimd_new(:,1)=(rprim(:,1)+rprim(:,2))*0.5_dp
 717            rprimd_new(:,2)=(rprim(:,2)-rprim(:,1))*0.5_dp
 718          case(3)
 719            rprimd_new(:,1)=(rprim(:,2)+rprim(:,3))*0.5_dp
 720            rprimd_new(:,2)=(rprim(:,1)+rprim(:,3))*0.5_dp
 721            rprimd_new(:,3)=(rprim(:,1)+rprim(:,2))*0.5_dp
 722          case(2)
 723            rprimd_new(:,1)=(-rprim(:,1)+rprim(:,2)+rprim(:,3))*0.5_dp
 724            rprimd_new(:,2)=( rprim(:,1)-rprim(:,2)+rprim(:,3))*0.5_dp
 725            rprimd_new(:,3)=( rprim(:,1)+rprim(:,2)-rprim(:,3))*0.5_dp
 726          case(7)
 727            rprimd_new(:,1)=( rprim(:,1)*2.0_dp+rprim(:,2)+rprim(:,3))/3.0_dp
 728            rprimd_new(:,2)=(-rprim(:,1)      +rprim(:,2)+rprim(:,3))/3.0_dp
 729            rprimd_new(:,3)=(-rprim(:,1)-rprim(:,2)*2.0_dp+rprim(:,3))/3.0_dp
 730          end select
 731          call symrelrot(nsym,rprimd,rprimd_new,symrel,tolsym)
 732 !        Produce xred in the new system of coordinates
 733          call xred2xcart(natrd,rprimd,xcart,xred)
 734          call xcart2xred(natrd,rprimd_new,xcart,xred)
 735 !        Produce tnons in the new system of coordinates
 736          ABI_MALLOC(tnons_cart,(3,nsym))
 737          call xred2xcart(nsym,rprimd,tnons_cart,tnons)
 738          call xcart2xred(nsym,rprimd_new,tnons_cart,tnons)
 739          ABI_FREE(tnons_cart)
 740 
 741          ! write(std_out,*)' after change of coordinates, nsym =',nsym
 742          ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
 743          ! do ii=1,nsym
 744          ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
 745          ! end do
 746 
 747          ! Prune the symmetry operations: suppress those with
 748          ! exactly the same point and magnetic part
 749          nsym_now=1
 750          do isym=2,nsym
 751            irreducible=1
 752            do jsym=1,nsym_now
 753              if(sum(abs(symrel(:,:,isym)-symrel(:,:,jsym)))==0 .and. symafm(isym)==symafm(jsym)) then
 754                irreducible=0
 755                exit
 756              end if
 757            end do
 758            if(irreducible==1)then
 759              nsym_now=nsym_now+1
 760              symrel(:,:,nsym_now)=symrel(:,:,isym)
 761              tnons(:,nsym_now)=tnons(:,isym)
 762              symafm(nsym_now)=symafm(isym)
 763            end if
 764          end do
 765          nsym=nsym_now
 766 !        Translate tnons in the ]-0.5,0.5] interval
 767          tnons(:,1:nsym)=tnons(:,1:nsym)-nint(tnons(:,1:nsym)-1.0d-8)
 768 
 769          ! DEBUG
 770          ! write(std_out,*)' after reduction, nsym =',nsym
 771          ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
 772          ! do ii=1,nsym
 773          ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
 774          ! end do
 775          ! ENDDEBUG
 776 
 777          ! Now that symrel, tnons and xred are expressed in the primitive
 778          ! axis system, update the geometric quantities
 779          rprimd(:,:)=rprimd_new(:,:)
 780          rprim(:,:)=rprimd_new(:,:)
 781          call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 782          call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
 783        end if
 784 
 785      end if
 786 
 787      if(natom/=natrd.and.multiplicity == 1)then
 788        ! Generate the full set of atoms from its knowledge in the irreducible part.
 789        call fillcell(chrgat,natom,natrd,nsym,nucdipmom,spinat,symafm,symrel,tnons,tolsym,typat,xred)
 790      end if
 791 
 792      ! Check whether the symmetry operations are consistent with the lattice vectors
 793      iexit=0
 794 
 795      call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tolsym)
 796 
 797    else
 798      ! spgroup==0 and nsym==0
 799 
 800      ! Here, spgroup==0 as well as nsym==0, so must generate
 801      ! the spatial group of symmetry. However, all the atom
 802      ! positions must be known, so the number
 803      ! of atoms to be read must equal the total number of atoms.
 804      if(natrd/=natom .and. multiplicity== 1)then
 805        write(msg, '(a,i0,a,a,i0,a,a,a,a,a,a,a,a,a)' )&
 806          'The number of atoms to be read (natrd)= ',natrd,ch10,&
 807          'differs from the total number of atoms (natom)= ',natom,ch10,&
 808          'while spgroup=0 and nsym=0.',&
 809          'This is not allowed, since the information needed to',ch10,&
 810          'generate the missing atomic coordinates is not available.',ch10,&
 811          'Action: modify your input file',ch10,&
 812          '(either natrd, or natom, or spgroup, or nsym)'
 813        ABI_ERROR(msg)
 814      endif
 815      if (multiplicity==1) typat(:)=typat_read(:)
 816 
 817      ! Find the symmetry operations: nsym, symafm, symrel and tnons.
 818      ! Use nptsym and ptsymrel, as determined by symlatt
 819 
 820      noncoll=0; if (nspden == 4) noncoll=1
 821 
 822      use_inversion=1
 823      if (dtset%usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then
 824        ABI_COMMENT("Removing inversion and improper rotations from initial space group because of PAW + SOC")
 825        use_inversion=0
 826      end if
 827 
 828      ! Get field in reduced coordinates (reduced e/d field)
 829 
 830      field_xred(:)=zero
 831      if (dtset%berryopt ==4) then
 832        do ii=1,3
 833          field_xred(ii)=dot_product(dtset%efield(:),gprimd(:,ii))
 834        end do
 835      else if (dtset%berryopt == 6 ) then
 836        do ii=1,3
 837          field_xred(ii)=dot_product(dtset%dfield(:),gprimd(:,ii))
 838          field_xred(ii)=field_xred(ii)+ dot_product(dtset%efield(:),gprimd(:,ii)) ! note: symmetry broken by D and E
 839        end do
 840      else if (dtset%berryopt == 14) then
 841        do ii=1,3
 842          field_xred(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii))
 843        end do
 844      else if (dtset%berryopt == 16) then
 845        do ii=1,3
 846          field_xred(ii)=dtset%red_dfield(ii)+dtset%red_efield(ii)  ! symmetry broken by reduced d and e
 847        end do
 848      else if (dtset%berryopt == 17) then
 849        do ii=1,3
 850          field_xred(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii))
 851          if(dtset%jfielddir(ii)==2) field_xred(ii)=dtset%red_dfield(ii)
 852        end do
 853      end if
 854 
 855      ! Loop on trials to generate better point symmetries by relying on a primitive cell instead (possibly) of a non-primitive one,
 856      ! This loop has been disactivated, because it is not clear that one can generate a more complete set of point symmetries 
 857      ! WITH INTEGER components of symrel from a primitive cell. One should allow non-integer components, but this would 
 858      ! be a large departure from the current implementation. Still, the detection of the existence of the primitive cell
 859      ! and the corresponding Bravais lattice is activated.
 860      do try_primitive=1,1
 861 
 862        call symfind(dtset%berryopt,field_xred,gprimd,jellslab,msym,natom,noncoll,nptsym,nsym,&
 863          nzchempot,dtset%prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
 864          chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr)
 865 
 866        !If the group closure is not obtained, which should be exceptional, try with a larger tolsym (three times larger)
 867        if(ierr/=0)then
 868          ABI_WARNING('Will try to obtain group closure by using a tripled tolsym.')
 869          call symfind(dtset%berryopt,field_xred,gprimd,jellslab,msym,natom,noncoll,nptsym,nsym,&
 870            nzchempot,dtset%prtvol,ptsymrel,spinat,symafm,symrel,tnons,three*tolsym,typat,use_inversion,xred,&
 871            chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr)
 872          ABI_CHECK(ierr==0,"Error in group closure")
 873          ABI_WARNING('Succeeded to obtain group closure by using a tripled tolsym.')
 874        endif
 875 
 876        ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize tnons for gliding or screw operations, 
 877        ! symmetrize the atomic positions and recompute the symmetry operations 
 878        if(tolsym>1.00001e-8)then
 879          call symmetrize_tnons(nsym,symrel,tnons,tolsym)
 880          ABI_MALLOC(indsym,(4,natom,nsym))
 881          ABI_MALLOC(symrec,(3,3,nsym))
 882          do isym=1,nsym
 883            call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
 884          end do
 885          call symatm(indsym,natom,nsym,symrec,tnons,tolsym,typat,xred)
 886          call symmetrize_xred(natom,nsym,symrel,tnons,xred,indsym=indsym)
 887          ABI_FREE(indsym)
 888          ABI_FREE(symrec)
 889 
 890          if(print_comment_tolsym==1)then
 891            write(msg,'(a,es12.3,18a)')&
 892 &           'The tolerance on symmetries =',tolsym,' is bigger than 1.0e-8.',ch10,&
 893 &           'In order to avoid spurious effects, the atomic coordinates have been',ch10,&
 894 &           'symmetrized before storing them in the dataset internal variable.',ch10,&
 895 &           'So, do not be surprised by the fact that your input variables (xcart, xred, ...)',ch10,&
 896 &           'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,&
 897 &           'This is not a problem per se.',ch10,&
 898 &           'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',&
 899 &           ' decrease tolsym to 1.0e-8 or lower,',ch10,&
 900 &           'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.',ch10,&
 901 &           'This message will only be printed once, even if there are other datasets where tolsym is bigger than 1.0e-8.'
 902            ABI_COMMENT(msg)
 903            print_comment_tolsym=0
 904          endif
 905 
 906          call symfind(dtset%berryopt,field_xred,gprimd,jellslab,msym,natom,noncoll,nptsym,nsym,&
 907            nzchempot,dtset%prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
 908            chrgat=chrgat,nucdipmom=nucdipmom)
 909 
 910          !Needs one more resymmetrization, for the tnons
 911          ABI_MALLOC(tnons_new,(3,nsym))
 912 
 913          call symmetrize_xred(natom,nsym,symrel,tnons,xred,&
 914 &          fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym)
 915          tnons(:,1:nsym)=tnons_new(:,:)
 916          ABI_FREE(tnons_new)
 917 
 918        end if
 919 
 920        chkprim_fake=-1 
 921        ABI_MALLOC(is_translation,(nsym))
 922        call chkprimit(chkprim_fake, multi, nsym, symafm, symrel, is_translation) 
 923 
 924        if(multi/=1)then ! The cell is not primitive, get the point symmetries from a primitive cell.
 925          ntranslat=multi
 926          ABI_MALLOC(translations,(3,ntranslat))
 927          itranslat=0
 928          do isym=1,nsym 
 929            if(is_translation(isym)==1)then
 930              itranslat=itranslat+1
 931              translations(:,itranslat)=tnons(:,isym)
 932            endif
 933          enddo
 934          ABI_FREE(is_translation)
 935          call reduce2primitive(ntranslat, rprimd, rprimd_primitive, tolsym, translations)
 936          ABI_FREE(translations)
 937          !Find the Bravais lattice of the primitive cell, and the point symmetries (however, in the primitive basis)
 938          call symlatt(bravais_reduced,msym,nptsym,ptsymrel,rprimd_primitive,tolsym)
 939          write(msg,'(2a,3(3es16.8,a),2(a,i4,a),3(a,3i4,a),a,i4)')&
 940 &          ' The cell is not primitive. One could obtain a primitive cell using the following primitive vectors (rprimd) :',ch10,&
 941 &          rprimd_primitive(1:3,1),ch10,&
 942 &          rprimd_primitive(1:3,2),ch10,&
 943 &          rprimd_primitive(1:3,3),ch10,&
 944 &          ' The Bravais lattice has iholohedry   =',bravais(1),ch10,&
 945 &          '                         center       =',bravais(2),ch10,&
 946 &          '                         bravais(3:5) =',bravais(3:5),ch10,&
 947 &          '                         bravais(6:8) =',bravais(6:8),ch10,&
 948 &          '                         bravais(9:11)=',bravais(9:11),ch10,&
 949 &          ' The number of point symmetries is nptsym=',nptsym           
 950          ABI_COMMENT(msg)
 951 
 952          !Convert the point symmetries to the non-primitive reduced coordinates
 953          call symrelrot(nsym, rprimd_primitive, rprimd, ptsymrel, tolsym, ierr)
 954          !Perhaps not all components of symrel are integers. This generates a return code, and precludes upgrading ptsymrel.
 955          if(ierr/=0)then
 956            write(msg,'(a)')&
 957 &           ' Not all components of symrel are integers in the primitive cell coordinate system.'
 958            ABI_COMMENT(msg)
 959            exit
 960          endif
 961        else ! The cell is primitive
 962          ABI_FREE(is_translation)
 963          exit
 964        endif
 965 
 966      enddo ! try_primitive
 967 
 968    end if
 969 
 970    ! Finalize the computation of coordinates: produce xcart
 971    call xred2xcart(natom,rprimd,xcart,xred)
 972 
 973  end if ! check of existence of an object
 974 
 975  ABI_FREE(ptsymrel)
 976  ABI_FREE(xcart_read)
 977  ABI_FREE(xcart)
 978  ABI_FREE(xred_read)
 979  ABI_FREE(typat_read)
 980 
 981  call geo%free()
 982 
 983  ! Correct the default nsym value, if a symmetry group has not been generated.
 984  if (nsym==0) nsym=1
 985 
 986 !--------------------------------------------------------------------------------------------------------
 987 
 988  icoulomb=0
 989  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'icoulomb',tread,'INT')
 990  if(tread==1)icoulomb=intarr(1)
 991 
 992  ! calculate the center of the atomic system such as to put the
 993  ! atoms in the middle of the simulation box for the free BC case.
 994  if (icoulomb == 1) then
 995    rcm(:)=zero
 996    do iatom=1,natom
 997      rcm(:)=rcm(:)+xred(:,iatom)
 998    end do
 999    rcm(:)=rcm(:)/real(natom,dp)-half
1000    do iatom=1,natom
1001      xred(:,iatom)=xred(:,iatom)-rcm(:)
1002    end do
1003    ! Also modify the tnons
1004    do isym=1,nsym
1005      tnons(:,isym)=matmul(symrel(:,:,isym),rcm(:))-rcm(:)+tnons(:,isym)
1006    end do
1007 
1008    ABI_WARNING('icoulomb is 1 --> the average center of coordinates has been translated to (0.5,0.5,0.5)')
1009  end if
1010 
1011  !========================================================================================================
1012  !
1013  ! At this stage, the cell parameters and atomic coordinates are known, as well as the symmetry operations
1014  ! There has been a preliminary analysis of the holohedry (not definitive, though ...)
1015  !
1016  !========================================================================================================
1017 
1018  ! Here, determine correctly the Bravais lattice and other space group or shubnikov group characteristics
1019  call symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym)
1020 
1021  ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize the rprimd. Keep xred fixed.
1022  if(tolsym>1.00001e-8)then
1023    ! Check whether the symmetry operations are consistent with the lattice vectors
1024    iexit=1
1025 
1026    call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tol8)
1027 
1028    if(iexit==-1)then
1029       write(msg,'(5a,es11.3,15a)')&
1030         'It is observed that the input primitive vectors are not accurate:',ch10,&
1031         'the lattice is not left invariant within 1.0e-8 when applying symmetry operations.',ch10,&
1032         'However, they are only slightly inaccurate, as inaccuracies are within the input tolsym=', tolsym,ch10,&
1033         'In order to avoid spurious effects, the primitive vectors have been',ch10,&
1034         'symmetrized before storing them in the dataset internal variable.',ch10,&
1035         'So, do not be surprised by the fact that your input variables (acell, rprim, xcart, xred, ...)',ch10,&
1036         'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,&
1037 &       'This is not a problem per se.',ch10,& 
1038 &       'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',&
1039 &       ' decrease tolsym to 1.0e-8 or lower.',ch10,&
1040         'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.'
1041      ABI_WARNING(msg)
1042 
1043      call symmetrize_rprimd(bravais,nsym,rprimd,symrel,tol8)
1044      call mkradim(acell,rprim,rprimd)
1045 
1046      !Needs one more resymmetrization, for the tnons
1047      ABI_MALLOC(tnons_new,(3,nsym))
1048      call symmetrize_xred(natom,nsym,symrel,tnons,xred,&
1049 &          fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym)
1050      tnons(:,1:nsym)=tnons_new(:,:)
1051      ABI_FREE(tnons_new)
1052 
1053    end if
1054 
1055  end if
1056 
1057  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1058  angdeg(1)=180.0_dp/pi * acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))
1059  angdeg(2)=180.0_dp/pi * acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))
1060  angdeg(3)=180.0_dp/pi * acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))
1061 !write(std_out,'(a,3f14.8)') ' ingeo: angdeg(1:3)=',angdeg(1:3)
1062 
1063 !--------------------------------------------------------------------------------------
1064 
1065  !Finally prune the set of symmetry in case non-symmorphic operations must be excluded
1066  if(symmorphi==0)then
1067    jsym=0
1068    do isym=1,nsym
1069      if(sum(tnons(:,isym)**2)<tol6)then
1070        jsym=jsym+1
1071        ! This symmetry operation is non-symmorphic, and can be kept
1072        if(isym/=jsym)then
1073          symrel(:,:,jsym)=symrel(:,:,isym)
1074          tnons(:,jsym)=tnons(:,isym)
1075          symafm(jsym)=symafm(isym)
1076        end if
1077      end if
1078    end do
1079    nsym=jsym
1080  end if
1081 
1082  !call symmultsg(nsym,symafm,symrel,tnons)
1083 
1084  ! 9) initialize the list of fixed atoms, and initial velocities -----------------
1085  ! Note: these inputs do not influence the previous generation of
1086  ! symmetry operations. This might be changed in the future
1087 
1088  ! idir=0 is for iatfix , idir=1 is for iatfixx,
1089  ! idir=2 is for iatfixy, idir=3 is for iatfixz
1090  iatfix(:,:)=0
1091 
1092  do idir=0,3
1093 
1094    if(idir==0)then
1095      call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfix',tread,'INT')
1096    else if(idir==1)then
1097      call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixx',tread,'INT')
1098    else if(idir==2)then
1099      call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixy',tread,'INT')
1100    else if(idir==3)then
1101      call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixz',tread,'INT')
1102    end if
1103 
1104    ! Use natfix also for natfixx,natfixy,natfixz
1105    natfix=0
1106    if(tread==1) natfix=intarr(1)
1107 
1108    ! Check the validity of natfix
1109    if (natfix<0 .or. natfix>natom) then
1110      write(msg, '(a,a,a,i0,a,i4,a,a,a)' )&
1111        'The input variables natfix, natfixx, natfixy and natfixz must be',ch10,&
1112        'between 0 and natom (= ',natom,'), while one of them is ',natfix,'.',ch10,&
1113        'Action: correct that occurence in your input file.'
1114      ABI_ERROR(msg)
1115    end if
1116 
1117    !Read iatfix
1118    if(idir==0)then
1119      call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfix',tread,'INT')
1120    else if(idir==1)then
1121      call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixx',tread,'INT')
1122    else if(idir==2)then
1123      call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixy',tread,'INT')
1124    else if(idir==3)then
1125      call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixz',tread,'INT')
1126    end if
1127 
1128    ! If some iatfix was read, natfix must vanish
1129    if (natfix==0 .and. tread==1)then
1130      write(msg, '(a,i1,5a)' )&
1131        'For direction ',idir,' the corresponding natfix is zero,',ch10,&
1132        'while iatfix specifies some atoms to be fixed.',ch10,&
1133        'Action: either specify a non-zero natfix(x,y,z) or suppress iatfix(x,y,z).'
1134      ABI_ERROR(msg)
1135    end if
1136 
1137    ! If natfix is non-zero, iatfix must be defined
1138    if (natfix>0 .and. tread==0)then
1139      write(msg, '(a,i1,3a,i0,3a)' )&
1140        'For direction ',idir,' no iatfix has been specified,',ch10,&
1141        'while natfix specifies that some atoms to be fixed, natfix= ',natfix,'.',ch10,&
1142        'Action: either set natfix(x,y,z) to zero or define iatfix(x,y,z).'
1143      ABI_ERROR(msg)
1144    end if
1145 
1146    if(tread==1)then
1147      do ii=1,natfix
1148        ! Checks the validity of the input iatfix
1149        if (intarr(ii)<1 .or. intarr(ii)>natom) then
1150          write(msg, '(a,a,a,i0,a,a,a)' )&
1151            'The input variables iatfix, iatfixx, iatfixy and iatfixz must be',ch10,&
1152            'between 1 and natom, while one of them is ',intarr(ii),'.',ch10,&
1153            'Action: correct that occurence in your input file.'
1154          ABI_ERROR(msg)
1155        end if
1156        ! Finally set the value of the internal iatfix array
1157        do iatom=1,natom
1158          if(intarr(ii)==iatom)then
1159            if(idir==0)iatfix(1:3,iatom)=1
1160            if(idir/=0)iatfix(idir,iatom)=1
1161          end if
1162        end do
1163      end do
1164    end if
1165 
1166  end do
1167 
1168  vel(:,:)=zero
1169  call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'vel',tread,'DPR')
1170  if(tread==1)vel(:,:)=reshape( dprarr(1:3*natom), [3, natom])
1171  call intagm_img(vel,iimage,jdtset,lenstr,nimage,3,natom,string,"vel",tread,'DPR')
1172 
1173  vel_cell(:,:)=zero
1174  call intagm(dprarr,intarr,jdtset,marr,3*3,string(1:lenstr),'vel_cell',tread,'DPR')
1175  if(tread==1)vel_cell(:,:)=reshape( dprarr(1:9), [3,3])
1176  call intagm_img(vel_cell,iimage,jdtset,lenstr,nimage,3,3,string,"vel_cell",tread,'DPR')
1177 
1178  ! mixalch
1179  if(ntypalch>0)then
1180    call intagm(dprarr,intarr,jdtset,marr,npspalch*ntypalch,string(1:lenstr),'mixalch',tread,'DPR')
1181    if(tread==1) mixalch(1:npspalch,1:ntypalch)= reshape(dprarr(1:npspalch*ntypalch), [npspalch, ntypalch])
1182    do itypat=1,ntypalch
1183      sumalch=sum(mixalch(1:npspalch,itypat))
1184      if(abs(sumalch-one)>tol10)then
1185        write(msg, '(a,i0,2a,f8.2,4a)' )&
1186          'For the alchemical atom number ',itypat,ch10,&
1187          'the sum of the pseudopotential coefficients is',sumalch,ch10,&
1188          'while it should be one.',ch10,&
1189          'Action: check the content of the input variable mixalch.'
1190        ABI_ERROR(msg)
1191      end if
1192    end do
1193    call intagm_img(mixalch,iimage,jdtset,lenstr,nimage,npspalch,ntypalch,string,"mixalch",tread,'DPR')
1194  end if
1195 
1196  ! amu (needs mixalch to be initialized ...)
1197  ! Find the default mass
1198  ABI_MALLOC(mass_psp,(npsp))
1199  do ipsp=1,npsp
1200    call atomdata_from_znucl(atom,znucl(ipsp))
1201    amu_default = atom%amu
1202    mass_psp(ipsp)=amu_default
1203  end do
1204  ! When the pseudo-atom is pure, simple copy
1205  ntyppure=ntypat-ntypalch
1206  if(ntyppure>0)then
1207    amu(1:ntyppure)=mass_psp(1:ntyppure)
1208  end if
1209  ! When the pseudo-atom is alchemical, must make mixing
1210  if(ntypalch>0)then
1211    do itypat=ntyppure+1,ntypat
1212      amu(itypat)=zero
1213      do ipsp=ntyppure+1,npsp
1214        amu(itypat)=amu(itypat)+mixalch(ipsp-ntyppure,itypat-ntyppure)*mass_psp(ipsp)
1215      end do
1216    end do
1217  end if
1218  ABI_FREE(mass_psp)
1219 
1220  call intagm(dprarr,intarr,jdtset,marr,ntypat,string(1:lenstr),'amu',tread,'DPR')
1221  if(tread==1)amu(:)=dprarr(1:ntypat)
1222  call intagm_img(amu,iimage,jdtset,lenstr,nimage,ntypat,string,"amu",tread,'DPR')
1223 
1224  ABI_FREE(intarr)
1225  ABI_FREE(dprarr)
1226 
1227 !DEBUG
1228 !write(std_out,'(a)')' m_ingeo%ingeo : exit '
1229 !call flush(std_out)
1230 !ENDDEBUG
1231 
1232 end subroutine ingeo

m_ingeo/ingeobld [ Functions ]

[ Top ] [ m_ingeo ] [ Functions ]

NAME

 ingeobld

FUNCTION

 The geometry builder.
 Start from the types and coordinates of the primitive atoms
 and produce the completed set of atoms, by using the definition
 of objects, then application of rotation, translation and repetition.

INPUTS

 iout=unit number of output file
 jdtset=number of the dataset looked for
 lenstr=actual length of the string
 natrd=number of atoms that have been read in the calling routine
 natom=number of atoms
 nobj=the number of objects
 string*(*)=character string containing all the input data. Initialized previously in instrng.
 typat_read(natrd)=type integer for each atom in the primitive set
 xcart_read(3,natrd)=cartesian coordinates of atoms (bohr), in the primitive set

OUTPUT

 typat(natom)=type integer for each atom in cell
 xcart(3,natom)=cartesian coordinates of atoms (bohr)

SOURCE

1262 subroutine ingeobld (iout,jdtset,lenstr,natrd,natom,nobj,string,typat,typat_read,xcart,xcart_read)
1263 
1264 !Arguments ------------------------------------
1265 !scalars
1266  integer,intent(in) :: iout,jdtset,lenstr,natom,natrd,nobj
1267  character(len=*),intent(in) :: string
1268 !arrays
1269  integer,intent(in) :: typat_read(natrd)
1270  integer,intent(out) :: typat(natom)
1271  real(dp),intent(in) :: xcart_read(3,natrd)
1272  real(dp),intent(out) :: xcart(3,natom)
1273 
1274 !Local variables-------------------------------
1275  character(len=*), parameter :: format01110 ="(1x,a6,1x,(t9,8i8) )"
1276  character(len=*), parameter :: format01160 ="(1x,a6,1x,1p,(t9,3g18.10)) "
1277 !scalars
1278  integer :: belonga,belongb,iatom,iatrd,ii,irep,irep1,irep2,irep3,ivac,marr
1279  integer :: natom_toberead,nread,objan,objbn,rotate,shift,tread,vacnum
1280  real(dp) :: angle,cosine,norm2per,norma,normb,normper,project,sine
1281  character(len=500) :: msg
1282 !arrays
1283  integer :: objarf(3),objbrf(3)
1284  integer,allocatable :: objaat(:),objbat(:),typat_full(:),vaclst(:)
1285  real(dp) :: axis2(3),axis3(3),axisa(3),axisb(3),objaax(6),objaro(4),objatr(12)
1286  real(dp) :: objbax(6),objbro(4),objbtr(12),parall(3),perpen(3),rotated(3)
1287  real(dp) :: vectora(3),vectorb(3)
1288  real(dp),allocatable :: xcart_full(:,:)
1289  integer,allocatable :: intarr(:)
1290  real(dp),allocatable :: dprarr(:)
1291 
1292 ! *************************************************************************
1293 
1294  marr=max(12,3*natom)
1295  ABI_MALLOC(intarr,(marr))
1296  ABI_MALLOC(dprarr,(marr))
1297 
1298 !1) Set up the number of vacancies.
1299 
1300 !This is the default
1301  vacnum=0
1302  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacnum',tread,'INT')
1303  if(tread==1) vacnum=intarr(1)
1304 
1305  if (vacnum>0)then
1306    ABI_MALLOC(vaclst,(vacnum))
1307 !  Read list of atoms to be suppressed to create vacancies
1308    call intagm(dprarr,intarr,jdtset,marr,vacnum,string(1:lenstr),'vaclst',tread,'INT')
1309    if(tread==1) vaclst(:)=intarr(1:vacnum)
1310    if(tread/=1)then
1311      write(msg, '(a,a,a,a,a)' )&
1312 &     'The array vaclst MUST be initialized in the input file',ch10,&
1313 &     'when vacnum is non-zero.',ch10,&
1314 &     'Action: initialize vaclst in your input file.'
1315      ABI_ERROR(msg)
1316    end if
1317  end if
1318 
1319  natom_toberead=natom+vacnum
1320 
1321 !2) Set up list and number of atoms in objects, and the --------------
1322 !operations to be performed on objects.
1323 
1324  write(msg,'(80a,a)')('=',ii=1,80),ch10
1325  call wrtout(std_out,msg)
1326  call wrtout(iout,msg)
1327 
1328  write(msg, '(a,a)' )'--ingeobld: echo values of variables connected to objects --------',ch10
1329  call wrtout(std_out,msg)
1330  call wrtout(iout,msg)
1331 
1332  if(vacnum>0)then
1333    write(iout,format01110) 'vacnum',vacnum
1334    write(std_out,format01110) 'vacnum',vacnum
1335    write(iout,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:)
1336    write(std_out,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:)
1337    write(iout, '(a)' ) ' '
1338    write(std_out,'(a)' ) ' '
1339  end if
1340 
1341  write(iout,format01110) 'nobj',nobj
1342  write(std_out,format01110) 'nobj',nobj
1343 
1344  if(nobj/=1 .and. nobj/=2)then
1345    write(msg, '(a,a,a,i8,a,a,a)' )&
1346 &   'The number of object (nobj) must be either 1 or 2,',ch10,&
1347 &   'while the input file has  nobj=',nobj,'.',ch10,&
1348 &   'Action: correct nobj in your input file.'
1349    ABI_ERROR(msg)
1350  end if
1351 
1352  if(nobj==1 .or. nobj==2)then
1353 
1354 !  Read the number of atoms of the object a
1355    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objan',tread,'INT')
1356    if(tread==1) objan=intarr(1)
1357 
1358    if(tread/=1)then
1359      write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1360 &     'The number of atoms in object a (objan) must be initialized',ch10,&
1361 &     'in the input file, when nobj=',nobj,'.',ch10,&
1362 &     'This is not the case.',ch10,&
1363 &     'Action: correct objan in your input file.'
1364      ABI_ERROR(msg)
1365    end if
1366 
1367    write(iout, '(a)' ) ' '
1368    write(std_out,'(a)' ) ' '
1369    write(iout,format01110) 'objan',objan
1370    write(std_out,format01110) 'objan',objan
1371 
1372    if(objan<=1 .or. objan>natom)then
1373      write(msg, '(a,a,a,a,a,i8,a,a,a)' )&
1374 &     'The number of atoms in object a (objan) must be larger than 0',ch10,&
1375 &     'and smaller than natom.',ch10,&
1376 &     'It is equal to ',objan,', an unacceptable value.',ch10,&
1377 &     'Action: correct objan in your input file.'
1378      ABI_ERROR(msg)
1379    end if
1380 
1381 !  Read list of atoms in object a
1382    call intagm(dprarr,intarr,jdtset,marr,objan,string(1:lenstr),'objaat',tread,'INT')
1383    ABI_MALLOC(objaat,(objan))
1384    if(tread==1) objaat(1:objan)=intarr(1:objan)
1385 
1386    if(tread/=1)then
1387      write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1388 &     'The list of atoms in object a (objaat) must be initialized',ch10,&
1389 &     'in the input file, when nobj=',nobj,'.',ch10,&
1390 &     'This is not the case.',ch10,&
1391 &     'Action: initialize objaat in your input file.'
1392      ABI_ERROR(msg)
1393    end if
1394 
1395    write(iout,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:)
1396    write(std_out,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:)
1397 
1398    do iatom=1,objan
1399      if(objaat(iatom)<1 .or. objaat(iatom)>natom)then
1400        write(msg, '(a,i8,a,a,i8,4a)' )&
1401 &       'The input value of objaat for atom number ',iatom,ch10,&
1402 &       'is equal to ',objaat(iatom),', an unacceptable value :',ch10,&
1403 &       'it should be between 1 and natom. ',&
1404 &       'Action: correct the array objaat in your input file.'
1405        ABI_ERROR(msg)
1406      end if
1407    end do
1408 
1409    if(objan>1)then
1410      do iatom=1,objan-1
1411        if( objaat(iatom)>=objaat(iatom+1) )then
1412          write(msg, '(a,i8,a,a,a,a,a,a)' )&
1413 &         'The input value of objaat for atom number ',iatom,ch10,&
1414 &         'is larger or equal to the one of the next atom,',ch10,&
1415 &         'while this list should be ordered, and an atom cannot be repeated.',ch10,&
1416 &         'Action: correct the array objaat in your input file.'
1417          ABI_ERROR(msg)
1418        end if
1419      end do
1420    end if
1421 
1422 !  Read repetition factors
1423    objarf(1:3)=1
1424    call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objarf',tread,'INT')
1425    if(tread==1) objarf(1:3)=intarr(1:3)
1426    write(iout,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:)
1427    write(std_out,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:)
1428 
1429    if(tread==1)then
1430      do irep=1,3
1431        if(objarf(irep)<1)then
1432          write(msg, '(a,a,a,3i8,a,a,a)' )&
1433 &         'The input values of objarf(1:3) must be positive,',ch10,&
1434 &         'while it is ',objarf(1:3),'.',ch10,&
1435 &         'Action: correct objarf in your input file.'
1436          ABI_ERROR(msg)
1437        end if
1438      end do
1439    end if
1440 
1441 !  Modify the number of atoms to be read
1442    natom_toberead=natom_toberead-objan*(objarf(1)*objarf(2)*objarf(3)-1)
1443 
1444 !  Read rotations angles and translations
1445    objaro(1:4)=0.0_dp
1446    objatr(1:12)=0.0_dp
1447    if (objarf(1)*objarf(2)*objarf(3) ==1) then
1448      nread=1
1449    else if (objarf(2)*objarf(3) ==1) then
1450      nread=2
1451    else if (objarf(3) ==1) then
1452      nread=3
1453    else
1454      nread=4
1455    end if
1456    call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objaro',tread,'DPR')
1457    if(tread==1) objaro(1:nread)=dprarr(1:nread)
1458 
1459    call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objatr',tread,'LEN')
1460 
1461    if(tread==1) objatr(1:3*nread)=dprarr(1:3*nread)
1462    write(iout,format01160) 'objaro',objaro(1:4)
1463    write(std_out,format01160) 'objaro',objaro(1:4)
1464    write(iout,format01160) 'objatr',objatr(1:12)
1465    write(std_out,format01160) 'objatr',objatr(1:12)
1466 !  If needed, read axes, but default to the x-axis to avoid errors later
1467    objaax(1:6)=0.0_dp ; objaax(4)=1.0_dp
1468 
1469    if(abs(objaro(1))+abs(objaro(2))+abs(objaro(3))+abs(objaro(4)) > 1.0d-10) then
1470      call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objaax',tread,'LEN')
1471      if(tread==1) objaax(1:6)=dprarr(1:6)
1472      if(tread/=1)then
1473        write(msg, '(a,a,a,a,a,a,a)' )&
1474 &       'The axis of object a (objaax) must be initialized',ch10,&
1475 &       'in the input file, when rotations (objaro) are present.',ch10,&
1476 &       'This is not the case.',ch10,&
1477 &       'Action: initialize objaax in your input file.'
1478        ABI_ERROR(msg)
1479      end if
1480      write(iout,format01160) 'objaax',objaax(1:6)
1481      write(std_out,format01160) 'objaax',objaax(1:6)
1482    end if
1483 
1484    axisa(1:3)=objaax(4:6)-objaax(1:3)
1485    norma=axisa(1)**2+axisa(2)**2+axisa(3)**2
1486 
1487    if(norma<1.0d-10)then
1488      write(msg, '(5a)' )&
1489 &     'The two points defined by the input array objaax are too',ch10,&
1490 &     'close to each other, and will not be used to define an axis.',ch10,&
1491 &     'Action: correct objaax in your input file.'
1492      ABI_ERROR(msg)
1493    end if
1494    axisa(1:3)=axisa(1:3)/sqrt(norma)
1495  end if !  End condition of existence of a first object
1496 
1497  if(nobj==2)then
1498 
1499 !  Read the number of atoms of the object b
1500    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objbn',tread,'INT')
1501    if(tread==1) objbn=intarr(1)
1502 
1503    if(tread/=1)then
1504      write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1505 &     'The number of atoms in object b (objbn) must be initialized',ch10,&
1506 &     'in the input file, when nobj=',nobj,'.',ch10,&
1507 &     'This is not the case.',ch10,&
1508 &     'Action: initialize objbn in your input file.'
1509      ABI_ERROR(msg)
1510    end if
1511 
1512    write(iout, '(a)' ) ' '
1513    write(std_out,'(a)' ) ' '
1514    write(iout,format01110) 'objbn',objbn
1515    write(std_out,format01110) 'objbn',objbn
1516 
1517    if(objbn<=1 .or. objbn>natom)then
1518      write(msg, '(a,a,a,a,a,i8,a,a,a)' )&
1519 &     'The number of atoms in object b (objbn) must be larger than 0',ch10,&
1520 &     'and smaller than natom.',ch10,&
1521 &     'It is equal to ',objbn,', an unacceptable value.',ch10,&
1522 &     'Action: correct objbn in your input file.'
1523      ABI_ERROR(msg)
1524    end if
1525 
1526 !  Read list of atoms in object b
1527    call intagm(dprarr,intarr,jdtset,marr,objbn,string(1:lenstr),'objbat',tread,'INT')
1528    ABI_MALLOC(objbat,(objbn))
1529 
1530    if(tread==1) objbat(1:objbn)=intarr(1:objbn)
1531    if(tread/=1)then
1532      write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1533 &     'The list of atoms in object b (objbat) must be initialized',ch10,&
1534 &     'in the input file, when nobj=',nobj,'.',ch10,&
1535 &     'This is not the case.',ch10,&
1536 &     'Action: initialize objbat in your input file.'
1537      ABI_ERROR(msg)
1538    end if
1539 
1540    write(iout,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:)
1541    write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:)
1542 
1543    do iatom=1,objbn
1544      if(objbat(iatom)<1 .or. objbat(iatom)>natom)then
1545        write(msg, '(a,i8,a,a,i8,a,a,a,a,a)' )&
1546 &       'The input value of objbat for atom number ',iatom,ch10,&
1547 &       'is equal to ',objbat(iatom),', an unacceptable value :',ch10,&
1548 &       'it should be between 1 and natom. ',ch10,&
1549 &       'Action: correct objbat in your input file.'
1550        ABI_ERROR(msg)
1551      end if
1552    end do
1553 
1554    if(objbn>1)then
1555      do iatom=1,objbn-1
1556        if( objbat(iatom)>=objbat(iatom+1) )then
1557          write(msg, '(a,i8,a,a,a,a,a,a)' )&
1558 &         'The input value of objbat for atom number ',iatom,ch10,&
1559 &         'is larger or equal to the one of the next atom,',ch10,&
1560 &         'while this list should be ordered, and an atom cannot be repeated.',ch10,&
1561 &         'Action: correct the array objbat in the input file.'
1562          ABI_ERROR(msg)
1563        end if
1564      end do
1565    end if
1566 
1567 !  Read repetition factors
1568    objbrf(1:3)=1
1569    call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objbrf',tread,'INT')
1570    if(tread==1) objbrf(1:3)=intarr(1:3)
1571    write(iout,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:)
1572    write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:)
1573 
1574    if(tread==1)then
1575      do irep=1,3
1576        if(objbrf(irep)<1)then
1577          write(msg, '(a,a,a,3i8,a,a,a)' )&
1578 &         'The input values of objbrf(1:3) must be positive,',ch10,&
1579 &         'while it is ',objbrf(1:3),'.',ch10,&
1580 &         'Action: correct objbrf in your input file.'
1581          ABI_ERROR(msg)
1582        end if
1583      end do
1584    end if
1585 
1586 !  Modify the number of atoms to be read
1587    natom_toberead=natom_toberead-objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1)
1588 !  Read rotations angles and translations
1589    objbro(1:4)=0.0_dp
1590    objbtr(1:12)=0.0_dp
1591    if (objbrf(1)*objbrf(2)*objbrf(3) ==1) then
1592      nread=1
1593    else if (objbrf(2)*objbrf(3) ==1) then
1594      nread=2
1595    else if (objbrf(3) ==1) then
1596      nread=3
1597    else
1598      nread=4
1599    end if
1600    call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objbro',tread,'DPR')
1601    if(tread==1) objbro(1:nread)=dprarr(1:nread)
1602 
1603    call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objbtr',tread,'LEN')
1604    if(tread==1) objbtr(1:3*nread)=dprarr(1:3*nread)
1605 
1606    write(iout,format01160) 'objbro',objbro(1:4)
1607    write(std_out,format01160) 'objbro',objbro(1:4)
1608    write(iout,format01160) 'objbtr',objbtr(1:12)
1609    write(std_out,format01160) 'objbtr',objbtr(1:12)
1610 
1611 !  If needed, read axes, but default to the x-axis to avoid errors later
1612    objbax(1:6)=0.0_dp ; objbax(4)=1.0_dp
1613    if(abs(objbro(1))+abs(objbro(2))+abs(objbro(3))+abs(objbro(4)) > 1.0d-10) then
1614      call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objbax',tread,'LEN')
1615      if(tread==1) objbax(1:6)=dprarr(1:6)
1616      if(tread/=1)then
1617        write(msg, '(a,a,a,a,a,a,a)' )&
1618 &       'The axis of object b (objbax) must be initialized',ch10,&
1619 &       'in the input file, when rotations (objbro) are present.',ch10,&
1620 &       'This is not the case.',ch10,&
1621 &       'Action: initialize objbax in your input file.'
1622        ABI_ERROR(msg)
1623      end if
1624      write(iout,format01160) 'objbax',objbax(1:6)
1625      write(std_out,format01160) 'objbax',objbax(1:6)
1626    end if
1627    axisb(1:3)=objbax(4:6)-objbax(1:3)
1628    normb=axisb(1)**2+axisb(2)**2+axisb(3)**2
1629    if(normb<1.0d-10)then
1630      write(msg, '(5a)' )&
1631 &     'The two points defined by the input array objbax are too',ch10,&
1632 &     'close to each other, and will not be used to define an axis.',ch10,&
1633 &     'Action: correct objbax in your input file.'
1634      ABI_ERROR(msg)
1635    end if
1636    axisb(1:3)=axisb(1:3)/sqrt(normb)
1637 
1638 !  Check whether both lists are disjoints. Use a very primitive algorithm.
1639    do iatom=1,objan
1640      do ii=1,objbn
1641        if(objaat(iatom)==objbat(ii))then
1642          write(msg, '(6a,i8,a,i8,3a)' )&
1643 &         'The objects a and b cannot have a common atom, but it is',ch10,&
1644 &         'found that the values of objaat and objbat ',&
1645 &         ' are identical, for their',ch10,&
1646 &         'atoms number ',iatom,' and ',ii,'.',ch10,&
1647 &         'Action: change objaat and/or objbat so that they have no common atom anymore.'
1648          ABI_ERROR(msg)
1649        end if
1650      end do
1651    end do
1652  end if !  End condition of existence of a second object
1653 
1654 !Check whether the number of atoms to be read obtained by relying
1655 !on natom, vacnum and the object definitions, or from natrd coincide
1656  if(natrd/=natom_toberead)then
1657    write(msg,'(11a,i0,a,i0,2a,i0,a)' )&
1658 &   ' ingeobld : ERROR -- ',ch10,&
1659 &   '  The number of atoms to be read (natrd) must be equal',ch10,&
1660 &   '  to the total number of atoms (natom), plus',ch10,&
1661 &   '  the number of vacancies (vacnum), minus',ch10,&
1662 &   '  the number of atoms added by the repetition of objects.',ch10,&
1663 &   '  This is not the case : natrd= ',natrd,', natom= ',natom,ch10,&
1664 &   ', vacnum= ',vacnum,';'
1665    call wrtout(std_out,msg)
1666 
1667    if(nobj==1 .or. nobj==2) then
1668      write(msg,'(a,i3,a,3i3,a,i5,a)' )&
1669 &     '   object a : objan=',objan,', objarf(1:3)=',objarf(1:3),&
1670 &     ' => adds ',objan*(objarf(1)*objarf(2)*objarf(3)-1),' atoms.'
1671      call wrtout(std_out,msg)
1672    end if
1673 
1674    if(nobj==2) then
1675      write(msg,'(a,i3,a,3i3,a,i5,a)' )&
1676 &     '   object b : objbn=',objbn,', objbrf(1:3)=',objbrf(1:3),&
1677 &     ' => adds ',objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1),' atoms.'
1678      call wrtout(std_out,msg)
1679    end if
1680 
1681    write(msg,'(3a)' )&
1682 &   '  Action : check the correspondence between natom+vacnum on one side,',ch10,&
1683 &   '           and natrd, objan, objbn, objarf and objbrf on the other side.'
1684    ABI_ERROR(msg)
1685  end if
1686 
1687 !6) Produce full set of atoms
1688 
1689 !Print the initial atom coordinates if the geometry builder is used
1690  write(iout, '(/,a)' )  ' Cartesian coordinates of the primitive atoms '
1691  write(std_out,'(/,a)' )' Cartesian coordinates of the primitive atoms '
1692  write(iout,format01160) '      ',xcart_read(:,:)
1693  write(std_out,format01160) '      ',xcart_read(:,:)
1694 
1695  ABI_MALLOC(typat_full,(natom+vacnum))
1696  ABI_MALLOC(xcart_full,(3,natom+vacnum))
1697 
1698 !Use the work array xcart_full to produce full set of atoms,
1699 !including those coming from repeated objects.
1700  iatom=1
1701  do iatrd=1,natrd
1702 
1703    belonga=0 ; belongb=0
1704    if(nobj==1 .or. nobj==2)then
1705 !    Determine whether the atom belongs to object a
1706      do ii=1,objan
1707        if(iatrd==objaat(ii))belonga=ii
1708      end do
1709    end if
1710    if(nobj==2)then
1711 !    Determine whether the atom belong to object b
1712      do ii=1,objbn
1713        if(iatrd==objbat(ii))belongb=ii
1714      end do
1715    end if
1716 
1717    !write(std_out,'(a,i5,a,i2,i2,a)' )' ingeobld : treating iatrd=',iatrd,', belong(a,b)=',belonga,belongb,'.'
1718 
1719 !  In case it does not belong to an object
1720    if(belonga==0 .and. belongb==0)then
1721      xcart_full(1:3,iatom)=xcart_read(1:3,iatrd)
1722      typat_full(iatom)=typat_read(iatrd)
1723      iatom=iatom+1
1724    else
1725 
1726 !    Repeat, rotate and translate this atom
1727      if(belonga/=0)then
1728 
1729 !      Treat object a
1730 !      Compute the relative coordinate of atom with respect to first point of axis
1731        vectora(1:3)=xcart_read(1:3,iatrd)-objaax(1:3)
1732 !      Project on axis
1733        project=vectora(1)*axisa(1)+vectora(2)*axisa(2)+vectora(3)*axisa(3)
1734 !      Get the parallel part
1735        parall(1:3)=project*axisa(1:3)
1736 !      Get the perpendicular part, to be rotated
1737        perpen(1:3)=vectora(1:3)-parall(1:3)
1738 !      Compute the norm of the perpendicular part
1739        norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2
1740 !      Initialisation to avoid warnings even if used behind if rotate == 1.
1741        normper = 0
1742 !      It the norm is too small, there is not need to rotate
1743        rotate=0
1744        if(norm2per>=1.0d-18)then
1745          rotate=1
1746          normper=sqrt(norm2per)
1747          axis2(1:3)=perpen(1:3)/normper
1748 !        Get the vector perpendicular to axisa and axisa2
1749          axis3(1)=axisa(2)*axis2(3)-axisa(3)*axis2(2)
1750          axis3(2)=axisa(3)*axis2(1)-axisa(1)*axis2(3)
1751          axis3(3)=axisa(1)*axis2(2)-axisa(2)*axis2(1)
1752        end if
1753 
1754 !      Here the repetition loop
1755        do irep3=1,objarf(3)
1756          do irep2=1,objarf(2)
1757            do irep1=1,objarf(1)
1758 !            Here the rotation
1759              if(rotate==1)then
1760 !              Compute the angle of rotation
1761                angle=objaro(1)+(irep1-1)*objaro(2) + &
1762 &               (irep2-1)*objaro(3)+(irep3-1)*objaro(4)
1763                cosine=cos(angle/180.0*pi)
1764                sine=sin(angle/180.0*pi)
1765                rotated(1:3)=objaax(1:3)+parall(1:3)+&
1766 &               normper*(cosine*axis2(1:3)+sine*axis3(1:3))
1767              else
1768                rotated(1:3)=vectora(1:3)
1769              end if
1770 !            Here the translation
1771              xcart_full(1:3,iatom)=rotated(1:3)+objatr(1:3)+&
1772 &             (irep1-1)*objatr(4:6)+(irep2-1)*objatr(7:9)+(irep3-1)*objatr(10:12)
1773              typat_full(iatom)=typat_read(iatrd)
1774              iatom=iatom+1
1775            end do
1776          end do
1777        end do ! End the repetition loop
1778 
1779      else
1780 !      If the atom belong to object b
1781 !      Compute the relative coordinate of atom with respect to first point of axis
1782        vectorb(1:3)=xcart_read(1:3,iatrd)-objbax(1:3)
1783 !      Project on axis
1784        project=vectorb(1)*axisb(1)+vectorb(2)*axisb(2)+vectorb(3)*axisb(3)
1785 !      Get the parallel part
1786        parall(1:3)=project*axisb(1:3)
1787 !      Get the perpendicular part, to be rotated
1788        perpen(1:3)=vectorb(1:3)-parall(1:3)
1789 !      Compute the norm of the perpendicular part
1790        norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2
1791 !      Initialisation to avoid warnings even if used behind if rotate == 1.
1792        normper = 0
1793 !      It the norm is too small, there is not need to rotate
1794        rotate=0
1795        if(norm2per>=1.0d-18)then
1796          rotate=1
1797          normper=sqrt(norm2per)
1798          axis2(1:3)=perpen(1:3)/normper
1799 !        Get the vector perpendicular to axisb and axis2
1800          axis3(1)=axisb(2)*axis2(3)-axisb(3)*axis2(2)
1801          axis3(2)=axisb(3)*axis2(1)-axisb(1)*axis2(3)
1802          axis3(3)=axisb(1)*axis2(2)-axisb(2)*axis2(1)
1803        end if
1804 !      Here the repetition loop
1805        do irep3=1,objbrf(3)
1806          do irep2=1,objbrf(2)
1807            do irep1=1,objbrf(1)
1808 !            Here the rotation
1809              if(rotate==1)then
1810 !              Compute the angle of rotation
1811                angle=objbro(1)+(irep1-1)*objbro(2) + &
1812 &               (irep2-1)*objbro(3)+ (irep3-1)*objbro(4)
1813                cosine=cos(angle/180.0*pi)
1814                sine=sin(angle/180.0*pi)
1815                rotated(1:3)=objbax(1:3)+parall(1:3)+&
1816 &               normper*(cosine*axis2(1:3)+sine*axis3(1:3))
1817              else
1818                rotated(1:3)=vectorb(1:3)
1819              end if
1820 !            Here the translation
1821              xcart_full(1:3,iatom)=rotated(1:3)+objbtr(1:3)+&
1822 &             (irep1-1)*objbtr(4:6)+(irep2-1)*objbtr(7:9)+(irep3-1)*objbtr(10:12)
1823              typat_full(iatom)=typat_read(iatrd)
1824              iatom=iatom+1
1825            end do
1826          end do
1827        end do ! End the repetition loop
1828      end if ! Condition of belonging to object b
1829    end if ! Condition of belonging to an object
1830  end do ! Loop on atoms
1831 
1832 !Create the vacancies here
1833  if(vacnum/=0)then
1834 !  First label the vacant atoms as belonging to typat 0
1835    do ivac=1,vacnum
1836      typat_full(vaclst(ivac))=0
1837    end do
1838 !  Then compact the arrays
1839    shift=0
1840    do iatom=1,natom
1841      if(typat_full(iatom+shift)==0) shift=shift+1
1842      if(shift/=0)then
1843        xcart_full(1:3,iatom)=xcart_full(1:3,iatom+shift)
1844        typat_full(iatom)=typat_full(iatom+shift)
1845      end if
1846    end do
1847  end if
1848 
1849 !Transfer the content of xcart_full and typat_full to the proper location
1850  xcart(:,1:natom)=xcart_full(:,1:natom)
1851  typat(1:natom)=typat_full(1:natom)
1852 
1853  ABI_FREE(typat_full)
1854  ABI_FREE(xcart_full)
1855  if(allocated(objaat)) then
1856    ABI_FREE(objaat)
1857  end if
1858  if(allocated(objbat)) then
1859    ABI_FREE(objbat)
1860  end if
1861 
1862  ABI_FREE(intarr)
1863  ABI_FREE(dprarr)
1864  if (vacnum>0)  then
1865    ABI_FREE(vaclst)
1866  end if
1867 
1868 end subroutine ingeobld

m_ingeo/invacuum [ Functions ]

[ Top ] [ m_ingeo ] [ Functions ]

NAME

 invacuum

FUNCTION

 Determine whether there is vacuum along some of the primitive directions in real space.

INPUTS

 jdtset=number of the dataset looked for
 lenstr=actual length of the string
 natom=number of atoms
 rprimd(3,3)=dimensional real space primitive translations (bohr)
 string*(*)=character string containing all the input data.
  Initialized previously in instrng.
 xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

 vacuum(3)= for each direction, 0 if no vacuum, 1 if vacuum

SOURCE

2053 subroutine invacuum(jdtset,lenstr,natom,rprimd,string,vacuum,xred)
2054 
2055 !Arguments ------------------------------------
2056 !scalars
2057  integer,intent(in) :: jdtset,lenstr,natom
2058  character(len=*),intent(in) :: string
2059 !arrays
2060  integer,intent(out) :: vacuum(3)
2061  real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
2062 
2063 !Local variables-------------------------------
2064 !scalars
2065  integer :: ia,ii,marr,tread
2066  real(dp) :: max_diff_xred,ucvol,vacwidth,vacxred
2067 !arrays
2068  integer,allocatable :: list(:)
2069  integer,allocatable :: intarr(:)
2070  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2071  real(dp),allocatable :: xred_sorted(:)
2072  real(dp),allocatable :: dprarr(:)
2073 
2074 ! *************************************************************************
2075 
2076 !Compute the maximum size of arrays intarr and dprarr
2077  marr=3
2078  ABI_MALLOC(intarr,(marr))
2079  ABI_MALLOC(dprarr,(marr))
2080 
2081 !Get metric quantities
2082  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2083 
2084 !Read vacwidth, or set the default
2085  vacwidth=10.0_dp
2086  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacwidth',tread,'LEN')
2087  if(tread==1) vacwidth=dprarr(1)
2088 
2089 !Read vacuum, or compute it using the atomic coordinates and vacwidth.
2090  vacuum(1:3)=0
2091  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'vacuum',tread,'INT')
2092 
2093  if(tread==1)then
2094    vacuum(1:3)=intarr(1:3)
2095  else
2096 !  For each direction, determine whether a vacuum space exists
2097    ABI_MALLOC(list,(natom))
2098    ABI_MALLOC(xred_sorted,(natom))
2099    do ii=1,3
2100 !    This is the minimum xred difference needed to have vacwidth
2101      vacxred=vacwidth*sqrt(sum(gprimd(:,ii)**2))
2102 !    Project the reduced coordinate in the [0.0_dp,1.0_dp[ interval
2103      xred_sorted(:)=mod(xred(ii,:),1.0_dp)
2104 !    list is dummy
2105      list(:)=0
2106 !    Sort xred_sorted
2107      call sort_dp(natom,xred_sorted,list,tol14)
2108      if(natom==1)then
2109        max_diff_xred=1.0_dp
2110      else
2111 !      Compute the difference between each pair of atom in the sorted order
2112        max_diff_xred=0.0_dp
2113        do ia=1,natom-1
2114          max_diff_xred=max(max_diff_xred,xred_sorted(ia+1)-xred_sorted(ia))
2115        end do
2116 !      Do not forget the image of the first atom in the next cell
2117        max_diff_xred=max(max_diff_xred,1.0_dp+xred_sorted(1)-xred_sorted(ia))
2118      end if
2119      if(vacxred<max_diff_xred+tol10)vacuum(ii)=1
2120    end do
2121    ABI_FREE(list)
2122    ABI_FREE(xred_sorted)
2123  end if
2124 
2125 !DEBUG
2126 !write(std_out,*)' invacuum : vacuum=',vacuum(1:3)
2127 !ENDDEBUG
2128 
2129  ABI_FREE(intarr)
2130  ABI_FREE(dprarr)
2131 
2132 end subroutine invacuum