TABLE OF CONTENTS


ABINIT/m_ingeo [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ingeo

FUNCTION

 Initialize geometry variables for the ABINIT code.

COPYRIGHT

  Copyright (C) 1998-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_ingeo
28 
29  use defs_basis
30  use defs_abitypes
31  use m_intagm_img
32  use m_abicore
33  use m_errors
34  use m_atomdata
35  use m_sort
36 
37  use m_symtk,      only : mati3inv, chkorthsy, symrelrot, mati3det, symmetrize_rprimd, symmetrize_xred, symatm
38  use m_spgbuilder, only : gensymspgr, gensymshub, gensymshub4
39  use m_symfind,    only : symfind, symanal, symlatt
40  use m_geometry,   only : mkradim, mkrdim, xcart2xred, xred2xcart, randomcellpos, metric
41  use m_parser,     only : intagm
42 
43  implicit none
44 
45  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

  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

PARENTS

      ingeo

CHILDREN

SOURCE

1874 subroutine fillcell(natom,natrd,nsym,nucdipmom,spinat,symafm,symrel,tnons,tolsym,typat,xred)
1875 
1876 
1877 !This section has been created automatically by the script Abilint (TD).
1878 !Do not modify the following lines by hand.
1879 #undef ABI_FUNC
1880 #define ABI_FUNC 'fillcell'
1881 !End of the abilint section
1882 
1883  implicit none
1884 
1885 !Arguments ------------------------------------
1886 !scalars
1887  integer,intent(in) :: natom,natrd,nsym
1888 !arrays
1889  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
1890  integer,intent(inout) :: typat(natom)
1891  real(dp),intent(in) :: tolsym
1892  real(dp),intent(in) :: tnons(3,nsym)
1893  real(dp),intent(inout) :: nucdipmom(3,natom),spinat(3,natom),xred(3,natom)
1894 
1895 !Local variables ------------------------------
1896 !scalars
1897  integer :: curat,flagch,flageq,ii,iij,jj,kk
1898  character(len=500) :: message
1899 !arrays
1900  integer :: bcktypat(nsym*natrd)
1901  real(dp) :: bckat(3),bcknucdipmom(3,nsym*natrd)
1902  real(dp) :: bckspinat(3,nsym*natrd),bckxred(3,nsym*natrd)
1903 
1904 ! *************************************************************************
1905 
1906 !DEBUG
1907 !write(std_out,*)' fillcell : enter with nsym, natrd= ',nsym,natrd
1908 !write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
1909 !do ii=1,nsym
1910 !write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
1911 !end do
1912 !write(std_out,*)' Describe the input atoms (index,typat,xred,spinat)'
1913 !do jj=1,natrd
1914 !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj)
1915 !end do
1916 !ENDDEBUG
1917 
1918  curat=0
1919 
1920 !Cycle over all the symmetry operations
1921  do ii=1,nsym
1922 
1923 !  Cycle over all the atoms in the assymetric unit cell
1924    do jj=1,natrd
1925 
1926 !    Symmetry operation application
1927      bckat(:)=matmul(symrel(:,:,ii),xred(:,jj))+tnons(:,ii)
1928 
1929 !    Normalization of the coordinates in [0,1)
1930      do iij=1,3
1931        do while (bckat(iij)<-tolsym)
1932          bckat(iij)=bckat(iij)+1.0d0
1933        end do
1934        do while (bckat(iij)>=1.0d0-tolsym)
1935          bckat(iij)=bckat(iij)-1.0d0
1936        end do
1937      end do
1938 
1939 !    Check for duplicate atoms
1940      flagch=0
1941      do kk=1,curat
1942        flageq=0
1943        if ( abs(bckxred(1,kk)-bckat(1))<tolsym  .and. &
1944 &       abs(bckxred(2,kk)-bckat(2))<tolsym  .and. &
1945 &       abs(bckxred(3,kk)-bckat(3))<tolsym       ) exit
1946        flagch=flagch+1
1947      end do
1948 
1949      if (flagch==curat) then
1950 !      Add the obtained atom to the bckxred list
1951        curat=curat+1
1952        bckxred(:,curat)=bckat
1953        bcktypat(curat)=typat(jj)
1954        bcknucdipmom(:,curat)=nucdipmom(:,jj)
1955        bckspinat(:,curat)=spinat(:,jj)*symafm(ii)
1956      end if
1957 
1958    end do
1959  end do
1960 
1961 !DEBUG
1962 !write(std_out,*)' fillcell : Proposed coordinates ='
1963 !do ii=1,curat
1964 !write(std_out,'(i4,3es16.6)' )ii,bckxred(:,ii)
1965 !end do
1966 !ENDDEBUG
1967 
1968  if (curat>natom) then
1969    write(message, '(a,i3,a,a,i7,a,a,a,a)' )&
1970 &   'The number of atoms obtained from symmetries, ',curat,ch10,&
1971 &   'is greater than the input number of atoms, natom=',natom,ch10,&
1972 &   'This is not allowed.',ch10,&
1973 &   'Action: modify natom or the symmetry data in the input file.'
1974    MSG_ERROR(message)
1975  end if
1976 
1977  if (curat<natom) then
1978    write(message, '(a,i3,a,a,i7,a,a,a,a)' )&
1979 &   'The number of atoms obtained from symmetries, ',curat,ch10,&
1980 &   'is lower than the input number of atoms, natom=',natom,ch10,&
1981 &   'This is not allowed.',ch10,&
1982 &   'Action: modify natom or the symmetry data in the input file.'
1983    MSG_ERROR(message)
1984  end if
1985 
1986 !Assignment of symmetry to xred
1987  xred(:,1:natom)=bckxred(:,1:natom)
1988  typat(1:natom)=bcktypat(1:natom)
1989  nucdipmom(1:3,1:natom)=bcknucdipmom(1:3,1:natom)
1990  spinat(1:3,1:natom)=bckspinat(1:3,1:natom)
1991 
1992 !DEBUG
1993 !write(std_out,*)' fillcell : exit with natom=',natom
1994 !write(std_out,*)' Describe the output atoms (index,typat,xred,spinat)'
1995 !do jj=1,natom
1996 !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj)
1997 !end do
1998 !ENDDEBUG
1999 
2000 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, used
  only if choice=1 or 3. Initialized previously in instrng.
 supercell_latt(3,3)=supercell lattice

OUTPUT

 acell(3)=length of primitive vectors
 amu(ntypat)=mass of each atomic type
 bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
 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!

PARENTS

      invars1

CHILDREN

      atomdata_from_znucl,chkorthsy,fillcell,gensymshub,gensymshub4
      gensymspgr,intagm_img,ingeobld,intagm,mati3inv,metric,mkradim,mkrdim
      randomcellpos,symanal,symatm,symfind,symlatt,symmetrize_rprimd
      symmetrize_xred,symrelrot,wrtout,xcart2xred,xred2xcart

SOURCE

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

PARENTS

      ingeo

CHILDREN

      intagm,wrtout

SOURCE

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

PARENTS

      invars1,invars2

CHILDREN

      intagm,metric,sort_dp

SOURCE

2031 subroutine invacuum(jdtset,lenstr,natom,rprimd,string,vacuum,xred)
2032 
2033 
2034 !This section has been created automatically by the script Abilint (TD).
2035 !Do not modify the following lines by hand.
2036 #undef ABI_FUNC
2037 #define ABI_FUNC 'invacuum'
2038 !End of the abilint section
2039 
2040  implicit none
2041 
2042 !Arguments ------------------------------------
2043 !scalars
2044  integer,intent(in) :: jdtset,lenstr,natom
2045  character(len=*),intent(in) :: string
2046 !arrays
2047  integer,intent(out) :: vacuum(3)
2048  real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
2049 
2050 !Local variables-------------------------------
2051 !scalars
2052  integer :: ia,ii,marr,tread
2053  real(dp) :: max_diff_xred,ucvol,vacwidth,vacxred
2054 !arrays
2055  integer,allocatable :: list(:)
2056  integer,allocatable :: intarr(:)
2057  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2058  real(dp),allocatable :: xred_sorted(:)
2059  real(dp),allocatable :: dprarr(:)
2060 
2061 ! *************************************************************************
2062 
2063 !Compute the maximum size of arrays intarr and dprarr
2064  marr=3
2065  ABI_ALLOCATE(intarr,(marr))
2066  ABI_ALLOCATE(dprarr,(marr))
2067 
2068 !Get metric quantities
2069  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2070 
2071 !Read vacwidth, or set the default
2072  vacwidth=10.0_dp
2073  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacwidth',tread,'LEN')
2074  if(tread==1) vacwidth=dprarr(1)
2075 
2076 !Read vacuum, or compute it using the atomic coordinates and vacwidth.
2077  vacuum(1:3)=0
2078  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'vacuum',tread,'INT')
2079 
2080  if(tread==1)then
2081    vacuum(1:3)=intarr(1:3)
2082  else
2083 !  For each direction, determine whether a vacuum space exists
2084    ABI_ALLOCATE(list,(natom))
2085    ABI_ALLOCATE(xred_sorted,(natom))
2086    do ii=1,3
2087 !    This is the minimum xred difference needed to have vacwidth
2088      vacxred=vacwidth*sqrt(sum(gprimd(:,ii)**2))
2089 !    Project the reduced coordinate in the [0.0_dp,1.0_dp[ interval
2090      xred_sorted(:)=mod(xred(ii,:),1.0_dp)
2091 !    list is dummy
2092      list(:)=0
2093 !    Sort xred_sorted
2094      call sort_dp(natom,xred_sorted,list,tol14)
2095      if(natom==1)then
2096        max_diff_xred=1.0_dp
2097      else
2098 !      Compute the difference between each pair of atom in the sorted order
2099        max_diff_xred=0.0_dp
2100        do ia=1,natom-1
2101          max_diff_xred=max(max_diff_xred,xred_sorted(ia+1)-xred_sorted(ia))
2102        end do
2103 !      Do not forget the image of the first atom in the next cell
2104        max_diff_xred=max(max_diff_xred,1.0_dp+xred_sorted(1)-xred_sorted(ia))
2105      end if
2106      if(vacxred<max_diff_xred+tol10)vacuum(ii)=1
2107    end do
2108    ABI_DEALLOCATE(list)
2109    ABI_DEALLOCATE(xred_sorted)
2110  end if
2111 
2112 !DEBUG
2113 !write(std_out,*)' invacuum : vacuum=',vacuum(1:3)
2114 !ENDDEBUG
2115 
2116  ABI_DEALLOCATE(intarr)
2117  ABI_DEALLOCATE(dprarr)
2118 
2119 end subroutine invacuum