TABLE OF CONTENTS
ABINIT/m_ingeo [ 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