TABLE OF CONTENTS
ABINIT/m_ingeo [ Modules ]
NAME
m_ingeo
FUNCTION
Initialize geometry variables for the ABINIT code.
COPYRIGHT
Copyright (C) 1998-2024 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, symfind_expert, 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
1891 subroutine fillcell(chrgat,natom,natrd,nsym,nucdipmom,spinat,symafm,symrel,tnons,tolsym,typat,xred) 1892 1893 !Arguments ------------------------------------ 1894 !scalars 1895 integer,intent(in) :: natom,natrd,nsym 1896 !arrays 1897 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 1898 integer,intent(inout) :: typat(natom) 1899 real(dp),intent(in) :: tolsym 1900 real(dp),intent(in) :: tnons(3,nsym) 1901 real(dp),intent(inout) :: chrgat(natom),nucdipmom(3,natom),spinat(3,natom),xred(3,natom) 1902 1903 !Local variables ------------------------------ 1904 !scalars 1905 integer :: curat,flagch,flageq,ii,iij,jj,kk 1906 character(len=500) :: msg 1907 !arrays 1908 integer :: bcktypat(nsym*natrd) 1909 real(dp) :: bckat(3),bcknucdipmom(3,nsym*natrd) 1910 real(dp) :: bckchrgat(nsym*natrd),bckspinat(3,nsym*natrd),bckxred(3,nsym*natrd) 1911 1912 ! ************************************************************************* 1913 1914 !DEBUG 1915 !write(std_out,*)' fillcell : enter with nsym, natrd= ',nsym,natrd 1916 !write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)' 1917 !do ii=1,nsym 1918 !write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii) 1919 !end do 1920 !write(std_out,*)' Describe the input atoms (index,typat,xred,spinat)' 1921 !do jj=1,natrd 1922 !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj) 1923 !end do 1924 !ENDDEBUG 1925 1926 curat=0 1927 1928 !Cycle over all the symmetry operations 1929 do ii=1,nsym 1930 1931 ! Cycle over all the atoms in the assymetric unit cell 1932 do jj=1,natrd 1933 1934 ! Symmetry operation application 1935 bckat(:)=matmul(symrel(:,:,ii),xred(:,jj))+tnons(:,ii) 1936 1937 ! Normalization of the coordinates in [0,1) 1938 do iij=1,3 1939 do while (bckat(iij)<-tolsym) 1940 bckat(iij)=bckat(iij)+1.0d0 1941 end do 1942 do while (bckat(iij)>=1.0d0-tolsym) 1943 bckat(iij)=bckat(iij)-1.0d0 1944 end do 1945 end do 1946 1947 ! Check for duplicate atoms 1948 flagch=0 1949 do kk=1,curat 1950 flageq=0 1951 if ( abs(bckxred(1,kk)-bckat(1))<tolsym .and. & 1952 & abs(bckxred(2,kk)-bckat(2))<tolsym .and. & 1953 & abs(bckxred(3,kk)-bckat(3))<tolsym ) exit 1954 flagch=flagch+1 1955 end do 1956 1957 if (flagch==curat) then 1958 ! Add the obtained atom to the bckxred list 1959 curat=curat+1 1960 bckxred(:,curat)=bckat 1961 bcktypat(curat)=typat(jj) 1962 bckchrgat(curat)=chrgat(jj) 1963 bcknucdipmom(:,curat)=nucdipmom(:,jj) 1964 bckspinat(:,curat)=spinat(:,jj)*symafm(ii) 1965 end if 1966 1967 end do 1968 end do 1969 1970 !DEBUG 1971 !write(std_out,*)' fillcell : Proposed coordinates =' 1972 !do ii=1,curat 1973 !write(std_out,'(i4,3es16.6)' )ii,bckxred(:,ii) 1974 !end do 1975 !ENDDEBUG 1976 1977 if (curat>natom) then 1978 write(msg, '(a,i3,a,a,i7,a,a,a,a)' )& 1979 & 'The number of atoms obtained from symmetries, ',curat,ch10,& 1980 & 'is greater 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 ABI_ERROR(msg) 1984 end if 1985 1986 if (curat<natom) then 1987 write(msg, '(a,i3,a,a,i7,a,a,a,a)' )& 1988 & 'The number of atoms obtained from symmetries, ',curat,ch10,& 1989 & 'is lower than the input number of atoms, natom=',natom,ch10,& 1990 & 'This is not allowed.',ch10,& 1991 & 'Action: modify natom or the symmetry data in the input file.' 1992 ABI_ERROR(msg) 1993 end if 1994 1995 !Assignment of symmetry to xred 1996 xred(:,1:natom)=bckxred(:,1:natom) 1997 typat(1:natom)=bcktypat(1:natom) 1998 chrgat(1:natom)=bckchrgat(1:natom) 1999 nucdipmom(1:3,1:natom)=bcknucdipmom(1:3,1:natom) 2000 spinat(1:3,1:natom)=bckspinat(1:3,1:natom) 2001 2002 !DEBUG 2003 !write(std_out,*)' fillcell : exit with natom=',natom 2004 !write(std_out,*)' Describe the output atoms (index,typat,xred,spinat)' 2005 !do jj=1,natom 2006 !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj) 2007 !end do 2008 !ENDDEBUG 2009 2010 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 field_red(3)=applied field direction in reduced coordinates 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
131 subroutine ingeo (acell,amu,bravais,chrgat,dtset,field_red,& 132 genafm,iatfix,icoulomb,iimage,iout,jdtset,jellslab,lenstr,mixalch,& 133 msym,natom,nimage,npsp,npspalch,nspden,nsppol,nsym,ntypalch,ntypat,& 134 nucdipmom,nzchempot,pawspnorb,& 135 ptgroupma,ratsph,rprim,slabzbeg,slabzend,spgroup,spinat,string,supercell_lattice,symafm,& 136 symmorphi,symrel,tnons,tolsym,typat,vel,vel_cell,xred,znucl,comm) 137 138 !Arguments ------------------------------------ 139 !scalars 140 integer,intent(in) :: iimage,iout,jdtset,lenstr,msym 141 integer,intent(in) :: nimage,npsp,npspalch,nspden,nsppol 142 integer,intent(in) :: ntypalch,ntypat,nzchempot,pawspnorb,comm 143 integer,intent(inout) :: natom,symmorphi 144 integer,intent(out) :: icoulomb,jellslab,ptgroupma,spgroup !vz_i 145 integer,intent(inout) :: nsym !vz_i 146 real(dp),intent(out) :: slabzbeg,slabzend,tolsym 147 character(len=*),intent(in) :: string 148 !arrays 149 integer,intent(in) :: supercell_lattice(3) 150 integer,intent(out) :: bravais(11),iatfix(3,natom) !vz_i 151 integer,intent(inout) :: symafm(msym) !vz_i 152 integer,intent(inout) :: symrel(3,3,msym) !vz_i 153 integer,intent(out) :: typat(natom) 154 real(dp),intent(inout) :: chrgat(natom) 155 real(dp),intent(inout) :: nucdipmom(3,natom) 156 real(dp),intent(in) :: ratsph(ntypat) 157 real(dp),intent(inout) :: spinat(3,natom) 158 real(dp),intent(out) :: acell(3),amu(ntypat),field_red(3),genafm(3),mixalch(npspalch,ntypalch) 159 real(dp),intent(inout) :: rprim(3,3),tnons(3,msym) !vz_i 160 real(dp),intent(out) :: vel(3,natom),vel_cell(3,3),xred(3,natom) 161 real(dp),intent(in) :: znucl(npsp) 162 type(dataset_type),intent(in) :: dtset 163 164 !Local variables------------------------------- 165 character(len=*), parameter :: format01110 ="(1x,a6,1x,(t9,8i8) )" 166 character(len=*), parameter :: format01160 ="(1x,a6,1x,1p,(t9,3g18.10)) " 167 !scalars 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 :: invar_z,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 :: intarr(:) 182 integer,allocatable :: is_translation(:) 183 integer,allocatable :: ptsymrel(:,:,:),typat_read(:) 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 if (dtset%berryopt ==4) then 257 do ii=1,3 258 field_red(ii)=dot_product(dtset%efield(:),gprimd(:,ii)) 259 end do 260 else if (dtset%berryopt == 6 ) then 261 do ii=1,3 262 field_red(ii)=dot_product(dtset%dfield(:),gprimd(:,ii)) 263 field_red(ii)=field_red(ii)+ dot_product(dtset%efield(:),gprimd(:,ii)) ! note: symmetry broken by D and E 264 end do 265 else if (dtset%berryopt == 14) then 266 do ii=1,3 267 field_red(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii)) 268 end do 269 else if (dtset%berryopt == 16) then 270 do ii=1,3 271 field_red(ii)=dtset%red_dfield(ii)+dtset%red_efield(ii) ! symmetry broken by reduced d and e 272 end do 273 else if (dtset%berryopt == 17) then 274 do ii=1,3 275 field_red(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii)) 276 if(dtset%jfielddir(ii)==2) field_red(ii)=dtset%red_dfield(ii) 277 end do 278 end if 279 280 !tolsym = tol8 281 !XG20200801 New default value for tolsym. This default value is also defined in m_invars1.F90 282 tolsym = tol5 283 !if (tread_geo /= 0 .and. geo%filetype == "poscar") then 284 ! tolsym = tol4 285 ! ABI_COMMENT("Reading structure from POSCAR --> default value of tolsym is set to 1e-4") 286 !end if 287 288 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'tolsym',tread,'DPR') 289 if(tread==1) tolsym=dprarr(1) 290 291 ! Find a tentative Bravais lattice and its point symmetries (might not use them) 292 ! Note that the Bravais lattice might not be the correct one yet (because the 293 ! actual atomic locations might lower the symmetry obtained from the lattice parameters only) 294 ABI_MALLOC(ptsymrel,(3,3,msym)) 295 call symlatt(bravais,dev_null,msym,nptsym,ptsymrel,rprimd,tolsym) 296 297 ! 3) Possibly, initialize a jellium slab 298 jellslab=0 299 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'jellslab',tread,'INT') 300 if(tread==1) jellslab=intarr(1) 301 302 slabzbeg=zero 303 slabzend=zero 304 if(jellslab/=0)then 305 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'slabzbeg',tread,'DPR') 306 if(tread==1) slabzbeg=dprarr(1) 307 308 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'slabzend',tread,'DPR') 309 if(tread==1) slabzend=dprarr(1) 310 end if 311 312 ! 4) Set up the number of atoms in the primitive set, to be read. 313 ! This is the default 314 natrd=natom 315 if(multiplicity > 1) natrd = natom_uc 316 317 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natrd',tnatrd,'INT') 318 if(tnatrd==1) natrd=intarr(1) 319 320 if(natrd<1 .or. natrd>natom)then 321 if(natrd>1 .and. multiplicity > 1) then 322 if(natrd < natom)then 323 write(msg, '(3a)' )& 324 'The number of atoms to be read (natrd) can not be used with supercell_latt.',ch10,& 325 'Action: Remove natrd or supercell_latt in your input file.' 326 ABI_ERROR(msg) 327 else 328 write(msg,'(3a,I0,a,I0,a,I0,2a)')& 329 'The input variable supercell_latt is present',ch10,& 330 'thus a supercell ',supercell_lattice(1),' ',supercell_lattice(2),& 331 ' ',supercell_lattice(3),' is generated',ch10 332 ABI_WARNING(msg) 333 end if 334 else 335 write(msg, '(3a,i0,a,i0,2a,a)' )& 336 'The number of atoms to be read (natrd) must be positive and not bigger than natom.',ch10,& 337 'This is not the case: natrd=',natrd,', natom=',natom,ch10,& 338 'Action: correct natrd or natom in your input file.' 339 ABI_ERROR(msg) 340 end if 341 end if 342 343 ! 5) Read the type and initial spin of each atom in the primitive set-------- 344 ABI_MALLOC(typat_read,(natrd)) 345 typat_read(1)=1 346 347 if (tread_geo == 0) then 348 call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'typat',tread,'INT') 349 350 ! If not read, try the XYZ data 351 if(tread==0) call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'_typat',tread,'INT') 352 if(tread==1) typat_read(1:natrd)=intarr(1:natrd) 353 354 else 355 typat_read = geo%typat 356 end if 357 358 do iatom=1,natrd 359 if(typat_read(iatom)<1 .or. typat_read(iatom)>ntypat )then 360 write(msg,'(a,i0,a,i0,a,a,a,i0,a,a,a)')& 361 'The input type of atom number ',iatom,' is equal to ',typat_read(iatom),',',ch10,& 362 'while it should be between 1 and ntypat= ',ntypat,'.',ch10,& 363 'Action: change either the variable typat or the variable ntypat.' 364 ABI_ERROR(msg) 365 end if 366 end do 367 368 ! 6) Read coordinates for each atom in the primitive set-------- 369 370 ABI_MALLOC(xcart_read,(3,natrd)) 371 ABI_MALLOC(xred_read,(3,natrd)) 372 373 random_atpos=0 374 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'random_atpos',txrandom,'INT') 375 if(txrandom==1) random_atpos=intarr(1) 376 if (random_atpos < 0 .or. random_atpos > 5) then 377 write(msg,'(3a)')& 378 'Random positions is a variable defined between 0 and 5. Error in the input file. ',ch10,& 379 'Action: define one of these in your input file.' 380 ABI_ERROR(msg) 381 end if 382 383 !if(nimage/=1 .and. iimage/=1)then 384 !FIXME: should this be called outside the above end if? 385 call randomcellpos(natom,npsp,ntypat,random_atpos,ratsph,rprim,rprimd_read,typat_read,& 386 xred_read(:,1:natrd),znucl,acell) 387 !This should not be printed if randomcellpos did nothing - it contains garbage. Spurious output anyway 388 !end if 389 390 if (tread_geo == 0) then 391 392 call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xred',txred,'DPR') 393 if (txred==1 .and. txrandom == 0) xred_read(:,1:natrd) = reshape(dprarr(1:3*natrd) , [3, natrd]) 394 call intagm_img(xred_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xred",txred,'DPR') 395 396 call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xcart',txcart,'LEN') 397 if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd]) 398 call intagm_img(xcart_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xcart",txcart,'LEN') 399 400 ! Might initialize xred from XYZ file 401 if (txred+txcart+txrandom==0) then 402 call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'_xred',txred,'DPR') 403 if (txred==1 .and. txrandom==0) xred_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd]) 404 405 call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'_xcart',txcart,'DPR') 406 if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd]) 407 end if 408 409 else 410 txcart = 0; txrandom = 0; txred = 1 411 xred_read = geo%xred 412 end if 413 414 if (txred + txcart + txrandom == 0) then 415 write(msg, '(3a)' )& 416 'Neither xred nor xcart are present in input file. ',ch10,& 417 'Action: define one of these in your input file.' 418 ABI_ERROR(msg) 419 end if 420 421 if (txred==1) write(msg, '(a)' ) ' xred is defined in input file' 422 if (txcart ==1) write(msg, '(a)' ) ' xcart is defined in input file (possibly in Angstrom)' 423 if (txrandom ==1) write(msg, '(a)' ) ' xred as random positions in the unit cell' 424 if (txrandom ==1) write(msg, '(a)' ) ' xcart are defined from a random distribution ' 425 call wrtout(std_out, msg) 426 427 if (txred + txcart + txrandom > 1)then 428 write(msg, '(3a)' )& 429 'Too many input channels for atomic positions are defined.',ch10,& 430 'Action: choose to define only one of these.' 431 ABI_ERROR(msg) 432 end if 433 434 if (txred==1 .or. txrandom /=0 ) then 435 call wrtout(std_out,' ingeo: takes atomic coordinates from input array xred ') 436 call xred2xcart(natrd,rprimd_read,xcart_read,xred_read) 437 else 438 call wrtout(std_out,' ingeo: takes atomic coordinates from input array xcart') 439 txcart=1 440 end if 441 442 !At this stage, the cartesian coordinates are known, for the atoms whose coordinates where read. 443 444 ! Here, allocate the variable that will contain the completed 445 ! sets of xcart, after the use of the geometry builder or the symmetry builder 446 ABI_MALLOC(xcart,(3,natom)) 447 448 !7) Eventually read the symmetries 449 !Take care of the symmetries 450 451 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nsym',tread,'INT') 452 if(tread==1) nsym=intarr(1) 453 454 ! Check that nsym is not negative 455 if (nsym<0) then 456 write(msg, '(a,i0,4a)' )& 457 'Input nsym must be positive or 0, but was ',nsym,ch10,& 458 'This is not allowed.',ch10,'Action: correct nsym in your input file.' 459 ABI_ERROR(msg) 460 end if 461 ! Check that nsym is not bigger than msym 462 if (nsym>msym) then 463 write(msg, '(2(a,i0),5a)')& 464 'Input nsym = ',nsym,' exceeds msym = ',msym,'.',ch10,& 465 'This is not allowed.',ch10,'Action: correct nsym in your input file.' 466 ABI_ERROR(msg) 467 end if 468 if (multiplicity>1) then 469 nsym = 1 470 ABI_WARNING('Input nsym is now set to one due to the supercell_latt input') 471 end if 472 473 ! Read symmorphi 474 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'symmorphi',tread,'INT') 475 if(tread==1) symmorphi=intarr(1) 476 477 ! Now, read the symmetry operations 478 if(nsym>0)then 479 call intagm(dprarr,intarr,jdtset,marr,9*nsym,string(1:lenstr),'symrel',tread,'INT') 480 if(nsym>1 .and. tread==0)then 481 write(msg,'(3a)')& 482 'When nsym>1, symrel must be defined in the input file.',ch10,& 483 'Action: either change nsym, or define symrel in your input file.' 484 ABI_ERROR(msg) 485 end if 486 if(tread==1) symrel(:,:,1:nsym)=reshape( intarr(1:9*nsym) , [3, 3, nsym]) 487 488 ! Take care of tnons 489 tnons(:,1:nsym)=zero 490 call intagm(dprarr,intarr,jdtset,marr,3*nsym,string(1:lenstr),'tnons',tread,'DPR') 491 if(tread==1) tnons(:,1:nsym)=reshape( dprarr(1:3*nsym), [3, nsym]) 492 493 if(symmorphi==0)then 494 do isym=1,nsym 495 if(sum(tnons(:,isym)**2)>tol6)then 496 write(msg, '(5a,i0,a,3f8.4,3a)' )& 497 'When symmorphi /= 1, the vectors of translation (tnons)',ch10,& 498 'a symmetry operation must vanish.',ch10,& 499 'However, for the symmetry operation number ',isym,', tnons =',tnons(:,isym),'.',ch10,& 500 'Action: either change your list of allowed symmetry operations, or use the symmetry finder (nsym=0).' 501 ABI_ERROR(msg) 502 end if 503 end do 504 end if 505 506 ! Take care of symafm 507 call intagm(dprarr,intarr,jdtset,marr,nsym,string(1:lenstr),'symafm',tread,'INT') 508 if(tread==1) symafm(1:nsym)=intarr(1:nsym) 509 end if 510 511 512 !8) Checks whether the geometry builder must be used, and call it if needed. 513 !Call the symmetry builder and analyzer if needed. 514 515 ! At this stage, nsym might still contain the default 0, msym contains the default dtset%maxnsym. 516 ! The cartesian coordinates of the atoms of the primitive set are contained in xcart_read. 517 518 nobj=0 519 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nobj',tread,'INT') 520 if(tread==1) nobj=intarr(1) 521 if(nobj /= 0 .and. multiplicity > 1)then 522 write(msg, '(3a)' )& 523 'nobj can not be used with supercell_latt.',ch10,& 524 'Action: Remove nobj or supercell_latt in your input file.' 525 ABI_ERROR(msg) 526 end if 527 528 !If there are objects, chkprim will not be used immediately 529 !But, if there are no objects, but a space group, it will be used directly. 530 !Need first to check the value of expert_user 531 expert_user=0 532 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'expert_user',tread,'INT') 533 if(tread==1) expert_user=intarr(1) 534 if(expert_user==0)then 535 chkprim=1 536 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'chkprim',tread,'INT') 537 if(tread==1) chkprim=intarr(1) 538 else 539 chkprim=0 540 endif 541 542 if(nobj/=0)then 543 544 ! chrgat is read for each atom, from 1 to natom 545 call intagm(dprarr,intarr,jdtset,marr,natom,string(1:lenstr),'chrgat',tread,'DPR') 546 if(tread==1) then 547 chrgat(1:natom) = dprarr(1:natom) 548 end if 549 550 ! Spinat is read for each atom, from 1 to natom 551 call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'spinat',tread,'DPR') 552 if(tread==1) then 553 spinat(1:3,1:natom) = reshape( dprarr(1:3*natom) , [3, natom]) 554 else if (nspden==4.or.(nspden==2.and.nsppol==1)) then 555 write(msg, '(5a)' )& 556 'When nspden=4 or (nspden==2 and nsppol==1), the input variable spinat must be',ch10,& 557 'defined in the input file, which is apparently not the case.',ch10,& 558 'Action: define spinat or use nspden=1 in your input file.' 559 ABI_ERROR(msg) 560 end if 561 562 ! nucdipmom is read for each atom, from 1 to natom 563 call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'nucdipmom',tread,'DPR') 564 if(tread==1) then 565 nucdipmom(1:3,1:natom) = reshape( dprarr(1:3*natom), [3, natom]) 566 end if 567 568 ! Will use the geometry builder 569 if(tnatrd/=1 .and. nobj/=0)then 570 write(msg, '(3a,i0,5a)' )& 571 'The number of atoms to be read (natrd) must be initialized',ch10,& 572 'in the input file, when nobj= ',nobj,'.',ch10,& 573 'This is not the case.',ch10,& 574 'Action: initialize natrd in your input file.' 575 ABI_ERROR(msg) 576 end if 577 578 if(jellslab/=0)then 579 write(msg, '(a,i0,3a)' )& 580 'A jellium slab cannot be used when nobj= ',nobj,'.',ch10,& 581 'Action: change one of the input variables jellslab or nobj in your input file.' 582 ABI_ERROR(msg) 583 end if 584 585 call ingeobld (iout,jdtset,lenstr,natrd,natom,nobj,string,typat,typat_read,xcart,xcart_read) 586 587 ! Finalize the computation of coordinates: produce xred. 588 call xcart2xred(natom,rprimd,xcart,xred) 589 590 else 591 ! nobj==0 592 593 ! chrgat is read for each irreducible atom, from 1 to natrd 594 call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'chrgat',tread,'DPR') 595 if(tread==1)chrgat(1:natrd) = dprarr(1:natrd) 596 597 ! Spinat is read for each irreducible atom, from 1 to natrd 598 call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'spinat',tread,'DPR') 599 if(tread==1)spinat(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natrd]) 600 601 ! nucdipmom is read for each irreducible atom, from 1 to natrd 602 call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'nucdipmom',tread,'DPR') 603 if(tread==1)nucdipmom(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natrd]) 604 605 ! Compute xred/typat and spinat for the supercell 606 if(multiplicity > 1)then 607 iatom_supercell = 0 608 do i1 = 1, supercell_lattice(1) 609 do i2 = 1, supercell_lattice(2) 610 do i3 = 1, supercell_lattice(3) 611 do iatom = 1, natom_uc 612 iatom_supercell = iatom_supercell + 1 613 xcart(:,iatom_supercell) = xcart_read(:,iatom) + matmul(rprimd_read,(/i1-1,i2-1,i3-1/)) 614 chrgat(iatom_supercell) = chrgat(iatom) 615 spinat(1:3,iatom_supercell) = spinat(1:3,iatom) 616 typat(iatom_supercell) = typat_read(iatom) 617 end do 618 end do 619 end do 620 end do 621 call xcart2xred(natom,rprimd,xcart,xred) 622 else 623 ! No supercell 624 call xcart2xred(natrd,rprimd,xcart_read,xred) 625 end if 626 627 spgroup=0 628 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgroup',tread,'INT') 629 if(tread==1) spgroup=intarr(1) 630 631 if(spgroup/=0 .or. nsym/=0)then 632 633 if(jellslab/=0 .and. nsym/=1 .and. spgroup/=1)then 634 write(msg, '(5a)' )& 635 'For the time being, a jellium slab can only be used',ch10,& 636 'either with the symmetry finder (nsym=0) or with the space group 1 (nsym=1)',ch10,& 637 'Action: change one of the input variables jellslab or nsym or spgroup in your input file.' 638 ABI_ERROR(msg) 639 end if 640 641 if(nzchempot/=0 .and. nsym/=1 .and. spgroup/=1)then 642 write(msg, '(5a)' )& 643 'For the time being, a spatially-varying chemical potential can only be used',ch10,& 644 'either with the symmetry finder (nsym=0) or with the space group 1 (nsym=1)',ch10,& 645 'Action: change one of the input variables nzchempot or nsym or spgroup in your input file.' 646 ABI_ERROR(msg) 647 end if 648 649 typat(1:natrd)=typat_read(1:natrd) 650 651 if(spgroup/=0 .and. nsym/=0)then 652 write(msg, '(a,i0,a,a,i0,a,a,a,a,a,a,a,a)' )& 653 'The spatial group number spgroup= ',spgroup,ch10,& 654 'is specified, as well as the number of symmetries nsym= ',nsym,ch10,& 655 'This is not allowed, as you can define the symmetries',ch10,& 656 'either using spgroup OR using nsym, but not both.',ch10,& 657 'Action: modify your input file',ch10,& 658 '(either set spgroup to 0, or nsym to 0)' 659 ABI_ERROR(msg) 660 end if 661 662 brvltt=0 663 664 if(spgroup/=0)then 665 666 ! Will generate the spatial group using spgroup 667 ! Assign default values 668 spgaxor=1 669 spgorig=1 670 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'brvltt',tread,'INT') 671 if(tread==1) brvltt=intarr(1) 672 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgaxor',tread,'INT') 673 if(tread==1) spgaxor=intarr(1) 674 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgorig',tread,'INT') 675 if(tread==1) spgorig=intarr(1) 676 677 ! Treat the case of magnetic groups 678 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgroupma',tspgroupma,'INT') 679 if(tspgroupma==1) spgroupma=intarr(1) 680 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'genafm',tgenafm,'DPR') 681 if(tgenafm==1) genafm(1:3)=dprarr(1:3) 682 if(tspgroupma/=0 .and. tgenafm/=0)then 683 write(msg, '(a,i0,a,a,3es9.2,a,a,a,a,a,a,a,a)' )& 684 'The spatial group number spgroupma= ',spgroupma,ch10,& 685 'is specified, as well as the antiferromagnetic generator genafm=',genafm(1:3),ch10,& 686 'This is not allowed, as you can define the magnetic space group',ch10,& 687 'either using spgroupma OR using genafm, but not both.',ch10,& 688 'Action: modify your input file',ch10,& 689 '(either define spgroupma or genafm)' 690 ABI_ERROR(msg) 691 end if 692 693 ! TODO: all the symmetry generation operations should be in one big routine 694 695 ! If spgroupma is defined, check whether it is consistent 696 ! with spgroup, determine the Shubnikov type, 697 ! and, for type IV, find the corresponding genafm 698 shubnikov=1 699 if(tspgroupma==1)then 700 call gensymshub(genafm,spgroup,spgroupma,shubnikov) 701 else if(tgenafm==1)then 702 shubnikov=4 703 end if 704 705 ! Generate the spatial group of symmetries in a conventional cell 706 ! In case of Shubnikov space group type IV, only generate the 707 ! Fedorov (non-magnetic) group. For Shubnikov type III space group, 708 ! the magnetic part is generated here. 709 bckbrvltt=brvltt 710 if(brvltt==-1)brvltt=0 711 call gensymspgr(brvltt,msym,nsym,shubnikov,spgaxor,spgorig,spgroup,spgroupma,symafm,symrel,tnons) 712 713 ! For shubnikov type IV groups, 714 ! double the space group, using the antiferromagnetic translation generator 715 if(shubnikov==4)then 716 call gensymshub4(genafm,msym,nsym,symafm,symrel,tnons) 717 end if 718 719 !write(std_out,*)' after gensymshub4, nsym =',nsym 720 !write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)' 721 !do ii=1,nsym 722 !write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii) 723 !end do 724 725 ! If brvltt was -1 at input, one should now change the conventional cell 726 ! to a primitive one, if brvltt/=1 727 if(bckbrvltt==-1 .and. brvltt/=1)then 728 ! Will work with rprim only 729 rprim(:,:)=rprimd(:,:) 730 rprimd_new(:,:)=rprimd(:,:) 731 acell(:)=1.0_dp 732 select case(brvltt) 733 case(5) 734 rprimd_new(:,2)=(rprim(:,2)+rprim(:,3))*0.5_dp 735 rprimd_new(:,3)=(rprim(:,3)-rprim(:,2))*0.5_dp 736 case(6) 737 rprimd_new(:,1)=(rprim(:,1)+rprim(:,3))*0.5_dp 738 rprimd_new(:,3)=(rprim(:,3)-rprim(:,1))*0.5_dp 739 case(4) 740 rprimd_new(:,1)=(rprim(:,1)+rprim(:,2))*0.5_dp 741 rprimd_new(:,2)=(rprim(:,2)-rprim(:,1))*0.5_dp 742 case(3) 743 rprimd_new(:,1)=(rprim(:,2)+rprim(:,3))*0.5_dp 744 rprimd_new(:,2)=(rprim(:,1)+rprim(:,3))*0.5_dp 745 rprimd_new(:,3)=(rprim(:,1)+rprim(:,2))*0.5_dp 746 case(2) 747 rprimd_new(:,1)=(-rprim(:,1)+rprim(:,2)+rprim(:,3))*0.5_dp 748 rprimd_new(:,2)=( rprim(:,1)-rprim(:,2)+rprim(:,3))*0.5_dp 749 rprimd_new(:,3)=( rprim(:,1)+rprim(:,2)-rprim(:,3))*0.5_dp 750 case(7) 751 rprimd_new(:,1)=( rprim(:,1)*2.0_dp+rprim(:,2)+rprim(:,3))/3.0_dp 752 rprimd_new(:,2)=(-rprim(:,1) +rprim(:,2)+rprim(:,3))/3.0_dp 753 rprimd_new(:,3)=(-rprim(:,1)-rprim(:,2)*2.0_dp+rprim(:,3))/3.0_dp 754 end select 755 call symrelrot(nsym,rprimd,rprimd_new,symrel,tolsym) 756 ! Produce xred in the new system of coordinates 757 call xred2xcart(natrd,rprimd,xcart,xred) 758 call xcart2xred(natrd,rprimd_new,xcart,xred) 759 ! Produce tnons in the new system of coordinates 760 ABI_MALLOC(tnons_cart,(3,nsym)) 761 call xred2xcart(nsym,rprimd,tnons_cart,tnons) 762 call xcart2xred(nsym,rprimd_new,tnons_cart,tnons) 763 ABI_FREE(tnons_cart) 764 765 ! write(std_out,*)' after change of coordinates, nsym =',nsym 766 ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)' 767 ! do ii=1,nsym 768 ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii) 769 ! end do 770 771 ! Prune the symmetry operations: suppress those with 772 ! exactly the same point and magnetic part 773 nsym_now=1 774 do isym=2,nsym 775 irreducible=1 776 do jsym=1,nsym_now 777 if(sum(abs(symrel(:,:,isym)-symrel(:,:,jsym)))==0 .and. symafm(isym)==symafm(jsym)) then 778 irreducible=0 779 exit 780 end if 781 end do 782 if(irreducible==1)then 783 nsym_now=nsym_now+1 784 symrel(:,:,nsym_now)=symrel(:,:,isym) 785 tnons(:,nsym_now)=tnons(:,isym) 786 symafm(nsym_now)=symafm(isym) 787 end if 788 end do 789 nsym=nsym_now 790 ! Translate tnons in the ]-0.5,0.5] interval 791 tnons(:,1:nsym)=tnons(:,1:nsym)-nint(tnons(:,1:nsym)-1.0d-8) 792 793 ! DEBUG 794 ! write(std_out,*)' after reduction, nsym =',nsym 795 ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)' 796 ! do ii=1,nsym 797 ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii) 798 ! end do 799 ! ENDDEBUG 800 801 ! Now that symrel, tnons and xred are expressed in the primitive 802 ! axis system, update the geometric quantities 803 rprimd(:,:)=rprimd_new(:,:) 804 rprim(:,:)=rprimd_new(:,:) 805 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 806 call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tolsym) 807 end if 808 809 end if 810 811 if(natom/=natrd.and.multiplicity == 1)then 812 ! Generate the full set of atoms from its knowledge in the irreducible part. 813 call fillcell(chrgat,natom,natrd,nsym,nucdipmom,spinat,symafm,symrel,tnons,tolsym,typat,xred) 814 end if 815 816 ! Check whether the symmetry operations are consistent with the lattice vectors 817 iexit=0 818 819 call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tolsym) 820 821 else 822 ! spgroup==0 and nsym==0 823 824 ! Here, spgroup==0 as well as nsym==0, so must generate 825 ! the spatial group of symmetry. However, all the atom 826 ! positions must be known, so the number 827 ! of atoms to be read must equal the total number of atoms. 828 if(natrd/=natom .and. multiplicity== 1)then 829 write(msg, '(a,i0,a,a,i0,a,a,a,a,a,a,a,a,a)' )& 830 'The number of atoms to be read (natrd)= ',natrd,ch10,& 831 'differs from the total number of atoms (natom)= ',natom,ch10,& 832 'while spgroup=0 and nsym=0.',& 833 'This is not allowed, since the information needed to',ch10,& 834 'generate the missing atomic coordinates is not available.',ch10,& 835 'Action: modify your input file',ch10,& 836 '(either natrd, or natom, or spgroup, or nsym)' 837 ABI_ERROR(msg) 838 endif 839 if (multiplicity==1) typat(:)=typat_read(:) 840 841 ! Find the symmetry operations: nsym, symafm, symrel and tnons. 842 ! Use nptsym and ptsymrel, as determined by symlatt 843 ! Will possibly correct xred and tnons. 844 845 noncoll=0; if (nspden == 4) noncoll=1 846 847 use_inversion=1 848 if (dtset%usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then 849 ABI_COMMENT("Removing inversion and improper rotations from initial space group because of PAW + SOC") 850 use_inversion=0 851 end if 852 853 ! Get field in reduced coordinates (reduced e/d field) 854 855 field_xred(:)=zero 856 if (dtset%berryopt ==4) then 857 do ii=1,3 858 field_xred(ii)=dot_product(dtset%efield(:),gprimd(:,ii)) 859 end do 860 else if (dtset%berryopt == 6 ) then 861 do ii=1,3 862 field_xred(ii)=dot_product(dtset%dfield(:),gprimd(:,ii)) 863 field_xred(ii)=field_xred(ii)+ dot_product(dtset%efield(:),gprimd(:,ii)) ! note: symmetry broken by D and E 864 end do 865 else if (dtset%berryopt == 14) then 866 do ii=1,3 867 field_xred(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii)) 868 end do 869 else if (dtset%berryopt == 16) then 870 do ii=1,3 871 field_xred(ii)=dtset%red_dfield(ii)+dtset%red_efield(ii) ! symmetry broken by reduced d and e 872 end do 873 else if (dtset%berryopt == 17) then 874 do ii=1,3 875 field_xred(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii)) 876 if(dtset%jfielddir(ii)==2) field_xred(ii)=dtset%red_dfield(ii) 877 end do 878 end if 879 880 ! Loop on trials to generate better point symmetries by relying on a primitive cell instead (possibly) of a non-primitive one, 881 ! This loop has been disactivated, because it is not clear that one can generate a more complete set of point symmetries 882 ! WITH INTEGER components of symrel from a primitive cell. One should allow non-integer components, but this would 883 ! be a large departure from the current implementation. Still, the detection of the existence of the primitive cell 884 ! and the corresponding Bravais lattice is activated. 885 do try_primitive=1,1 886 887 invar_z=0 ; if(jellslab/=0 .or. nzchempot/=0)invar_z=2 888 call symfind_expert(gprimd,msym,natom,nptsym,nspden,nsym,& 889 pawspnorb,dtset%prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,dtset%usepaw,xred,& 890 chrgat=chrgat,nucdipmom=nucdipmom,invardir_red=dtset%field_red,invar_z=invar_z) 891 892 chkprim_fake=-1 893 ABI_MALLOC(is_translation,(nsym)) 894 call chkprimit(chkprim_fake, multi, nsym, symafm, symrel, is_translation) 895 896 if(multi/=1)then ! The cell is not primitive, get the point symmetries from a primitive cell. 897 ntranslat=multi 898 ABI_MALLOC(translations,(3,ntranslat)) 899 itranslat=0 900 do isym=1,nsym 901 if(is_translation(isym)==1)then 902 itranslat=itranslat+1 903 translations(:,itranslat)=tnons(:,isym) 904 endif 905 enddo 906 ABI_FREE(is_translation) 907 call reduce2primitive(ntranslat, rprimd, rprimd_primitive, tolsym, translations) 908 ABI_FREE(translations) 909 !Find the Bravais lattice of the primitive cell, and the point symmetries (however, in the primitive basis) 910 call symlatt(bravais_reduced,dev_null,msym,nptsym,ptsymrel,rprimd_primitive,tolsym) 911 write(msg,'(2a,3(3es16.8,a),2(a,i4,a),3(a,3i4,a),a,i4)')& 912 & ' The cell is not primitive. One could obtain a primitive cell using the following primitive vectors (rprimd) :',ch10,& 913 & rprimd_primitive(1:3,1),ch10,& 914 & rprimd_primitive(1:3,2),ch10,& 915 & rprimd_primitive(1:3,3),ch10,& 916 & ' This Bravais lattice has iholohedry =',bravais(1),ch10,& 917 & ' center =',bravais(2),ch10,& 918 & ' bravais(3:5) =',bravais(3:5),ch10,& 919 & ' bravais(6:8) =',bravais(6:8),ch10,& 920 & ' bravais(9:11)=',bravais(9:11),ch10,& 921 & ' The number of point symmetries would be nptsym=',nptsym 922 ABI_COMMENT(msg) 923 924 !Convert the point symmetries to the non-primitive reduced coordinates 925 call symrelrot(nsym, rprimd_primitive, rprimd, ptsymrel, tolsym, ierr) 926 !Perhaps not all components of symrel are integers. This generates a return code, and precludes upgrading ptsymrel. 927 if(ierr/=0)then 928 write(msg,'(a)')& 929 & ' Not all components of symrel are integers in the primitive cell coordinate system.' 930 ABI_COMMENT(msg) 931 exit 932 endif 933 else ! The cell is primitive 934 ABI_FREE(is_translation) 935 exit 936 endif 937 938 enddo ! try_primitive 939 940 end if ! spgroup==0 and nsym==0 941 942 ! Finalize the computation of coordinates: produce xcart 943 call xred2xcart(natom,rprimd,xcart,xred) 944 945 end if ! check of existence of an object 946 947 ABI_FREE(ptsymrel) 948 ABI_FREE(xcart_read) 949 ABI_FREE(xcart) 950 ABI_FREE(xred_read) 951 ABI_FREE(typat_read) 952 953 call geo%free() 954 955 ! Correct the default nsym value, if a symmetry group has not been generated. 956 if (nsym==0) nsym=1 957 958 !-------------------------------------------------------------------------------------------------------- 959 960 icoulomb=0 961 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'icoulomb',tread,'INT') 962 if(tread==1)icoulomb=intarr(1) 963 964 ! calculate the center of the atomic system such as to put the 965 ! atoms in the middle of the simulation box for the free BC case. 966 if (icoulomb == 1) then 967 rcm(:)=zero 968 do iatom=1,natom 969 rcm(:)=rcm(:)+xred(:,iatom) 970 end do 971 rcm(:)=rcm(:)/real(natom,dp)-half 972 do iatom=1,natom 973 xred(:,iatom)=xred(:,iatom)-rcm(:) 974 end do 975 ! Also modify the tnons 976 do isym=1,nsym 977 tnons(:,isym)=matmul(symrel(:,:,isym),rcm(:))-rcm(:)+tnons(:,isym) 978 end do 979 980 ABI_WARNING('icoulomb is 1 --> the average center of coordinates has been translated to (0.5,0.5,0.5)') 981 end if 982 983 !======================================================================================================== 984 ! 985 ! At this stage, the cell parameters and atomic coordinates are known, as well as the symmetry operations 986 ! There has been a preliminary analysis of the holohedry (not definitive, though ...) 987 ! 988 !======================================================================================================== 989 990 !DEBUG 991 !write(std_out,'(a)')' m_ingeo%ingeo : before symanal ' 992 !call flush(std_out) 993 !ENDDEBUG 994 995 ! Here, determine correctly the Bravais lattice and other space group or shubnikov group characteristics 996 call symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym) 997 998 !DEBUG 999 !write(std_out,'(a)')' m_ingeo%ingeo : after symanal ' 1000 !call flush(std_out) 1001 !ENDDEBUG 1002 1003 ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize the rprimd. Keep xred fixed. 1004 if(tolsym>1.00001e-8)then 1005 ! Check whether the symmetry operations are consistent with the lattice vectors 1006 iexit=1 1007 1008 call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tol8) 1009 1010 if(iexit==-1)then 1011 write(msg,'(5a,es11.3,15a)')& 1012 'It is observed that the input primitive vectors are not accurate:',ch10,& 1013 'the lattice is not left invariant within 1.0e-8 when applying symmetry operations.',ch10,& 1014 'However, they are only slightly inaccurate, as inaccuracies are within the input tolsym=', tolsym,ch10,& 1015 'In order to avoid spurious effects, the primitive vectors have been',ch10,& 1016 'symmetrized before storing them in the dataset internal variable.',ch10,& 1017 'So, do not be surprised by the fact that your input variables (acell, rprim, xcart, xred, ...)',ch10,& 1018 'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,& 1019 & 'This is not a problem per se.',ch10,& 1020 & 'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',& 1021 & ' decrease tolsym to 1.0e-8 or lower.',ch10,& 1022 'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.' 1023 ABI_WARNING(msg) 1024 1025 call symmetrize_rprimd(bravais,nsym,rprimd,symrel,tol8) 1026 call mkradim(acell,rprim,rprimd) 1027 1028 !Needs one more resymmetrization, for the tnons 1029 ABI_MALLOC(tnons_new,(3,nsym)) 1030 call symmetrize_xred(natom,nsym,symrel,tnons,xred,& 1031 & fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym) 1032 tnons(:,1:nsym)=tnons_new(:,:) 1033 ABI_FREE(tnons_new) 1034 1035 end if 1036 1037 end if 1038 1039 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1040 angdeg(1)=180.0_dp/pi * acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3))) 1041 angdeg(2)=180.0_dp/pi * acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3))) 1042 angdeg(3)=180.0_dp/pi * acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2))) 1043 !write(std_out,'(a,3f14.8)') ' ingeo: angdeg(1:3)=',angdeg(1:3) 1044 1045 !-------------------------------------------------------------------------------------- 1046 1047 !Finally prune the set of symmetry in case non-symmorphic operations must be excluded 1048 if(symmorphi==0)then 1049 jsym=0 1050 do isym=1,nsym 1051 if(sum(tnons(:,isym)**2)<tol6)then 1052 jsym=jsym+1 1053 ! This symmetry operation is non-symmorphic, and can be kept 1054 if(isym/=jsym)then 1055 symrel(:,:,jsym)=symrel(:,:,isym) 1056 tnons(:,jsym)=tnons(:,isym) 1057 symafm(jsym)=symafm(isym) 1058 end if 1059 end if 1060 end do 1061 nsym=jsym 1062 end if 1063 1064 !call symmultsg(nsym,symafm,symrel,tnons) 1065 1066 ! 9) initialize the list of fixed atoms, and initial velocities ----------------- 1067 ! Note: these inputs do not influence the previous generation of 1068 ! symmetry operations. This might be changed in the future 1069 1070 ! idir=0 is for iatfix , idir=1 is for iatfixx, 1071 ! idir=2 is for iatfixy, idir=3 is for iatfixz 1072 iatfix(:,:)=0 1073 1074 do idir=0,3 1075 1076 if(idir==0)then 1077 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfix',tread,'INT') 1078 else if(idir==1)then 1079 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixx',tread,'INT') 1080 else if(idir==2)then 1081 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixy',tread,'INT') 1082 else if(idir==3)then 1083 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixz',tread,'INT') 1084 end if 1085 1086 ! Use natfix also for natfixx,natfixy,natfixz 1087 natfix=0 1088 if(tread==1) natfix=intarr(1) 1089 1090 ! Check the validity of natfix 1091 if (natfix<0 .or. natfix>natom) then 1092 write(msg, '(a,a,a,i0,a,i4,a,a,a)' )& 1093 'The input variables natfix, natfixx, natfixy and natfixz must be',ch10,& 1094 'between 0 and natom (= ',natom,'), while one of them is ',natfix,'.',ch10,& 1095 'Action: correct that occurence in your input file.' 1096 ABI_ERROR(msg) 1097 end if 1098 1099 !Read iatfix 1100 if(idir==0)then 1101 call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfix',tread,'INT') 1102 else if(idir==1)then 1103 call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixx',tread,'INT') 1104 else if(idir==2)then 1105 call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixy',tread,'INT') 1106 else if(idir==3)then 1107 call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixz',tread,'INT') 1108 end if 1109 1110 ! If some iatfix was read, natfix must vanish 1111 if (natfix==0 .and. tread==1)then 1112 write(msg, '(a,i1,5a)' )& 1113 'For direction ',idir,' the corresponding natfix is zero,',ch10,& 1114 'while iatfix specifies some atoms to be fixed.',ch10,& 1115 'Action: either specify a non-zero natfix(x,y,z) or suppress iatfix(x,y,z).' 1116 ABI_ERROR(msg) 1117 end if 1118 1119 ! If natfix is non-zero, iatfix must be defined 1120 if (natfix>0 .and. tread==0)then 1121 write(msg, '(a,i1,3a,i0,3a)' )& 1122 'For direction ',idir,' no iatfix has been specified,',ch10,& 1123 'while natfix specifies that some atoms to be fixed, natfix= ',natfix,'.',ch10,& 1124 'Action: either set natfix(x,y,z) to zero or define iatfix(x,y,z).' 1125 ABI_ERROR(msg) 1126 end if 1127 1128 if(tread==1)then 1129 do ii=1,natfix 1130 ! Checks the validity of the input iatfix 1131 if (intarr(ii)<1 .or. intarr(ii)>natom) then 1132 write(msg, '(a,a,a,i0,a,a,a)' )& 1133 'The input variables iatfix, iatfixx, iatfixy and iatfixz must be',ch10,& 1134 'between 1 and natom, while one of them is ',intarr(ii),'.',ch10,& 1135 'Action: correct that occurence in your input file.' 1136 ABI_ERROR(msg) 1137 end if 1138 ! Finally set the value of the internal iatfix array 1139 do iatom=1,natom 1140 if(intarr(ii)==iatom)then 1141 if(idir==0)iatfix(1:3,iatom)=1 1142 if(idir/=0)iatfix(idir,iatom)=1 1143 end if 1144 end do 1145 end do 1146 end if 1147 1148 end do 1149 1150 vel(:,:)=zero 1151 call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'vel',tread,'DPR') 1152 if(tread==1)vel(:,:)=reshape( dprarr(1:3*natom), [3, natom]) 1153 call intagm_img(vel,iimage,jdtset,lenstr,nimage,3,natom,string,"vel",tread,'DPR') 1154 1155 vel_cell(:,:)=zero 1156 call intagm(dprarr,intarr,jdtset,marr,3*3,string(1:lenstr),'vel_cell',tread,'DPR') 1157 if(tread==1)vel_cell(:,:)=reshape( dprarr(1:9), [3,3]) 1158 call intagm_img(vel_cell,iimage,jdtset,lenstr,nimage,3,3,string,"vel_cell",tread,'DPR') 1159 1160 ! mixalch 1161 if(ntypalch>0)then 1162 call intagm(dprarr,intarr,jdtset,marr,npspalch*ntypalch,string(1:lenstr),'mixalch',tread,'DPR') 1163 if(tread==1) mixalch(1:npspalch,1:ntypalch)= reshape(dprarr(1:npspalch*ntypalch), [npspalch, ntypalch]) 1164 do itypat=1,ntypalch 1165 sumalch=sum(mixalch(1:npspalch,itypat)) 1166 if(abs(sumalch-one)>tol10)then 1167 write(msg, '(a,i0,2a,f8.2,4a)' )& 1168 'For the alchemical atom number ',itypat,ch10,& 1169 'the sum of the pseudopotential coefficients is',sumalch,ch10,& 1170 'while it should be one.',ch10,& 1171 'Action: check the content of the input variable mixalch.' 1172 ABI_ERROR(msg) 1173 end if 1174 end do 1175 call intagm_img(mixalch,iimage,jdtset,lenstr,nimage,npspalch,ntypalch,string,"mixalch",tread,'DPR') 1176 end if 1177 1178 ! amu (needs mixalch to be initialized ...) 1179 ! Find the default mass 1180 ABI_MALLOC(mass_psp,(npsp)) 1181 do ipsp=1,npsp 1182 call atomdata_from_znucl(atom,znucl(ipsp)) 1183 amu_default = atom%amu 1184 mass_psp(ipsp)=amu_default 1185 end do 1186 ! When the pseudo-atom is pure, simple copy 1187 ntyppure=ntypat-ntypalch 1188 if(ntyppure>0)then 1189 amu(1:ntyppure)=mass_psp(1:ntyppure) 1190 end if 1191 ! When the pseudo-atom is alchemical, must make mixing 1192 if(ntypalch>0)then 1193 do itypat=ntyppure+1,ntypat 1194 amu(itypat)=zero 1195 do ipsp=ntyppure+1,npsp 1196 amu(itypat)=amu(itypat)+mixalch(ipsp-ntyppure,itypat-ntyppure)*mass_psp(ipsp) 1197 end do 1198 end do 1199 end if 1200 ABI_FREE(mass_psp) 1201 1202 call intagm(dprarr,intarr,jdtset,marr,ntypat,string(1:lenstr),'amu',tread,'DPR') 1203 if(tread==1)amu(:)=dprarr(1:ntypat) 1204 call intagm_img(amu,iimage,jdtset,lenstr,nimage,ntypat,string,"amu",tread,'DPR') 1205 1206 ABI_FREE(intarr) 1207 ABI_FREE(dprarr) 1208 1209 !DEBUG 1210 !write(std_out,'(a)')' m_ingeo%ingeo : exit ' 1211 !call flush(std_out) 1212 !ENDDEBUG 1213 1214 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
1244 subroutine ingeobld (iout,jdtset,lenstr,natrd,natom,nobj,string,typat,typat_read,xcart,xcart_read) 1245 1246 !Arguments ------------------------------------ 1247 !scalars 1248 integer,intent(in) :: iout,jdtset,lenstr,natom,natrd,nobj 1249 character(len=*),intent(in) :: string 1250 !arrays 1251 integer,intent(in) :: typat_read(natrd) 1252 integer,intent(out) :: typat(natom) 1253 real(dp),intent(in) :: xcart_read(3,natrd) 1254 real(dp),intent(out) :: xcart(3,natom) 1255 1256 !Local variables------------------------------- 1257 character(len=*), parameter :: format01110 ="(1x,a6,1x,(t9,8i8) )" 1258 character(len=*), parameter :: format01160 ="(1x,a6,1x,1p,(t9,3g18.10)) " 1259 !scalars 1260 integer :: belonga,belongb,iatom,iatrd,ii,irep,irep1,irep2,irep3,ivac,marr 1261 integer :: natom_toberead,nread,objan,objbn,rotate,shift,tread,vacnum 1262 real(dp) :: angle,cosine,norm2per,norma,normb,normper,project,sine 1263 character(len=500) :: msg 1264 !arrays 1265 integer :: objarf(3),objbrf(3) 1266 integer,allocatable :: objaat(:),objbat(:),typat_full(:),vaclst(:) 1267 real(dp) :: axis2(3),axis3(3),axisa(3),axisb(3),objaax(6),objaro(4),objatr(12) 1268 real(dp) :: objbax(6),objbro(4),objbtr(12),parall(3),perpen(3),rotated(3) 1269 real(dp) :: vectora(3),vectorb(3) 1270 real(dp),allocatable :: xcart_full(:,:) 1271 integer,allocatable :: intarr(:) 1272 real(dp),allocatable :: dprarr(:) 1273 1274 ! ************************************************************************* 1275 1276 marr=max(12,3*natom) 1277 ABI_MALLOC(intarr,(marr)) 1278 ABI_MALLOC(dprarr,(marr)) 1279 1280 !1) Set up the number of vacancies. 1281 1282 !This is the default 1283 vacnum=0 1284 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacnum',tread,'INT') 1285 if(tread==1) vacnum=intarr(1) 1286 1287 if (vacnum>0)then 1288 ABI_MALLOC(vaclst,(vacnum)) 1289 ! Read list of atoms to be suppressed to create vacancies 1290 call intagm(dprarr,intarr,jdtset,marr,vacnum,string(1:lenstr),'vaclst',tread,'INT') 1291 if(tread==1) vaclst(:)=intarr(1:vacnum) 1292 if(tread/=1)then 1293 write(msg, '(a,a,a,a,a)' )& 1294 & 'The array vaclst MUST be initialized in the input file',ch10,& 1295 & 'when vacnum is non-zero.',ch10,& 1296 & 'Action: initialize vaclst in your input file.' 1297 ABI_ERROR(msg) 1298 end if 1299 end if 1300 1301 natom_toberead=natom+vacnum 1302 1303 !2) Set up list and number of atoms in objects, and the -------------- 1304 !operations to be performed on objects. 1305 1306 write(msg,'(80a,a)')('=',ii=1,80),ch10 1307 call wrtout(std_out,msg) 1308 call wrtout(iout,msg) 1309 1310 write(msg, '(a,a)' )'--ingeobld: echo values of variables connected to objects --------',ch10 1311 call wrtout(std_out,msg) 1312 call wrtout(iout,msg) 1313 1314 if(vacnum>0)then 1315 write(iout,format01110) 'vacnum',vacnum 1316 write(std_out,format01110) 'vacnum',vacnum 1317 write(iout,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:) 1318 write(std_out,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:) 1319 write(iout, '(a)' ) ' ' 1320 write(std_out,'(a)' ) ' ' 1321 end if 1322 1323 write(iout,format01110) 'nobj',nobj 1324 write(std_out,format01110) 'nobj',nobj 1325 1326 if(nobj/=1 .and. nobj/=2)then 1327 write(msg, '(a,a,a,i8,a,a,a)' )& 1328 & 'The number of object (nobj) must be either 1 or 2,',ch10,& 1329 & 'while the input file has nobj=',nobj,'.',ch10,& 1330 & 'Action: correct nobj in your input file.' 1331 ABI_ERROR(msg) 1332 end if 1333 1334 if(nobj==1 .or. nobj==2)then 1335 1336 ! Read the number of atoms of the object a 1337 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objan',tread,'INT') 1338 if(tread==1) objan=intarr(1) 1339 1340 if(tread/=1)then 1341 write(msg, '(a,a,a,i8,a,a,a,a,a)' )& 1342 & 'The number of atoms in object a (objan) must be initialized',ch10,& 1343 & 'in the input file, when nobj=',nobj,'.',ch10,& 1344 & 'This is not the case.',ch10,& 1345 & 'Action: correct objan in your input file.' 1346 ABI_ERROR(msg) 1347 end if 1348 1349 write(iout, '(a)' ) ' ' 1350 write(std_out,'(a)' ) ' ' 1351 write(iout,format01110) 'objan',objan 1352 write(std_out,format01110) 'objan',objan 1353 1354 if(objan<=1 .or. objan>natom)then 1355 write(msg, '(a,a,a,a,a,i8,a,a,a)' )& 1356 & 'The number of atoms in object a (objan) must be larger than 0',ch10,& 1357 & 'and smaller than natom.',ch10,& 1358 & 'It is equal to ',objan,', an unacceptable value.',ch10,& 1359 & 'Action: correct objan in your input file.' 1360 ABI_ERROR(msg) 1361 end if 1362 1363 ! Read list of atoms in object a 1364 call intagm(dprarr,intarr,jdtset,marr,objan,string(1:lenstr),'objaat',tread,'INT') 1365 ABI_MALLOC(objaat,(objan)) 1366 if(tread==1) objaat(1:objan)=intarr(1:objan) 1367 1368 if(tread/=1)then 1369 write(msg, '(a,a,a,i8,a,a,a,a,a)' )& 1370 & 'The list of atoms in object a (objaat) must be initialized',ch10,& 1371 & 'in the input file, when nobj=',nobj,'.',ch10,& 1372 & 'This is not the case.',ch10,& 1373 & 'Action: initialize objaat in your input file.' 1374 ABI_ERROR(msg) 1375 end if 1376 1377 write(iout,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:) 1378 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:) 1379 1380 do iatom=1,objan 1381 if(objaat(iatom)<1 .or. objaat(iatom)>natom)then 1382 write(msg, '(a,i8,a,a,i8,4a)' )& 1383 & 'The input value of objaat for atom number ',iatom,ch10,& 1384 & 'is equal to ',objaat(iatom),', an unacceptable value :',ch10,& 1385 & 'it should be between 1 and natom. ',& 1386 & 'Action: correct the array objaat in your input file.' 1387 ABI_ERROR(msg) 1388 end if 1389 end do 1390 1391 if(objan>1)then 1392 do iatom=1,objan-1 1393 if( objaat(iatom)>=objaat(iatom+1) )then 1394 write(msg, '(a,i8,a,a,a,a,a,a)' )& 1395 & 'The input value of objaat for atom number ',iatom,ch10,& 1396 & 'is larger or equal to the one of the next atom,',ch10,& 1397 & 'while this list should be ordered, and an atom cannot be repeated.',ch10,& 1398 & 'Action: correct the array objaat in your input file.' 1399 ABI_ERROR(msg) 1400 end if 1401 end do 1402 end if 1403 1404 ! Read repetition factors 1405 objarf(1:3)=1 1406 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objarf',tread,'INT') 1407 if(tread==1) objarf(1:3)=intarr(1:3) 1408 write(iout,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:) 1409 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:) 1410 1411 if(tread==1)then 1412 do irep=1,3 1413 if(objarf(irep)<1)then 1414 write(msg, '(a,a,a,3i8,a,a,a)' )& 1415 & 'The input values of objarf(1:3) must be positive,',ch10,& 1416 & 'while it is ',objarf(1:3),'.',ch10,& 1417 & 'Action: correct objarf in your input file.' 1418 ABI_ERROR(msg) 1419 end if 1420 end do 1421 end if 1422 1423 ! Modify the number of atoms to be read 1424 natom_toberead=natom_toberead-objan*(objarf(1)*objarf(2)*objarf(3)-1) 1425 1426 ! Read rotations angles and translations 1427 objaro(1:4)=0.0_dp 1428 objatr(1:12)=0.0_dp 1429 if (objarf(1)*objarf(2)*objarf(3) ==1) then 1430 nread=1 1431 else if (objarf(2)*objarf(3) ==1) then 1432 nread=2 1433 else if (objarf(3) ==1) then 1434 nread=3 1435 else 1436 nread=4 1437 end if 1438 call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objaro',tread,'DPR') 1439 if(tread==1) objaro(1:nread)=dprarr(1:nread) 1440 1441 call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objatr',tread,'LEN') 1442 1443 if(tread==1) objatr(1:3*nread)=dprarr(1:3*nread) 1444 write(iout,format01160) 'objaro',objaro(1:4) 1445 write(std_out,format01160) 'objaro',objaro(1:4) 1446 write(iout,format01160) 'objatr',objatr(1:12) 1447 write(std_out,format01160) 'objatr',objatr(1:12) 1448 ! If needed, read axes, but default to the x-axis to avoid errors later 1449 objaax(1:6)=0.0_dp ; objaax(4)=1.0_dp 1450 1451 if(abs(objaro(1))+abs(objaro(2))+abs(objaro(3))+abs(objaro(4)) > 1.0d-10) then 1452 call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objaax',tread,'LEN') 1453 if(tread==1) objaax(1:6)=dprarr(1:6) 1454 if(tread/=1)then 1455 write(msg, '(a,a,a,a,a,a,a)' )& 1456 & 'The axis of object a (objaax) must be initialized',ch10,& 1457 & 'in the input file, when rotations (objaro) are present.',ch10,& 1458 & 'This is not the case.',ch10,& 1459 & 'Action: initialize objaax in your input file.' 1460 ABI_ERROR(msg) 1461 end if 1462 write(iout,format01160) 'objaax',objaax(1:6) 1463 write(std_out,format01160) 'objaax',objaax(1:6) 1464 end if 1465 1466 axisa(1:3)=objaax(4:6)-objaax(1:3) 1467 norma=axisa(1)**2+axisa(2)**2+axisa(3)**2 1468 1469 if(norma<1.0d-10)then 1470 write(msg, '(5a)' )& 1471 & 'The two points defined by the input array objaax are too',ch10,& 1472 & 'close to each other, and will not be used to define an axis.',ch10,& 1473 & 'Action: correct objaax in your input file.' 1474 ABI_ERROR(msg) 1475 end if 1476 axisa(1:3)=axisa(1:3)/sqrt(norma) 1477 end if ! End condition of existence of a first object 1478 1479 if(nobj==2)then 1480 1481 ! Read the number of atoms of the object b 1482 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objbn',tread,'INT') 1483 if(tread==1) objbn=intarr(1) 1484 1485 if(tread/=1)then 1486 write(msg, '(a,a,a,i8,a,a,a,a,a)' )& 1487 & 'The number of atoms in object b (objbn) must be initialized',ch10,& 1488 & 'in the input file, when nobj=',nobj,'.',ch10,& 1489 & 'This is not the case.',ch10,& 1490 & 'Action: initialize objbn in your input file.' 1491 ABI_ERROR(msg) 1492 end if 1493 1494 write(iout, '(a)' ) ' ' 1495 write(std_out,'(a)' ) ' ' 1496 write(iout,format01110) 'objbn',objbn 1497 write(std_out,format01110) 'objbn',objbn 1498 1499 if(objbn<=1 .or. objbn>natom)then 1500 write(msg, '(a,a,a,a,a,i8,a,a,a)' )& 1501 & 'The number of atoms in object b (objbn) must be larger than 0',ch10,& 1502 & 'and smaller than natom.',ch10,& 1503 & 'It is equal to ',objbn,', an unacceptable value.',ch10,& 1504 & 'Action: correct objbn in your input file.' 1505 ABI_ERROR(msg) 1506 end if 1507 1508 ! Read list of atoms in object b 1509 call intagm(dprarr,intarr,jdtset,marr,objbn,string(1:lenstr),'objbat',tread,'INT') 1510 ABI_MALLOC(objbat,(objbn)) 1511 1512 if(tread==1) objbat(1:objbn)=intarr(1:objbn) 1513 if(tread/=1)then 1514 write(msg, '(a,a,a,i8,a,a,a,a,a)' )& 1515 & 'The list of atoms in object b (objbat) must be initialized',ch10,& 1516 & 'in the input file, when nobj=',nobj,'.',ch10,& 1517 & 'This is not the case.',ch10,& 1518 & 'Action: initialize objbat in your input file.' 1519 ABI_ERROR(msg) 1520 end if 1521 1522 write(iout,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:) 1523 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:) 1524 1525 do iatom=1,objbn 1526 if(objbat(iatom)<1 .or. objbat(iatom)>natom)then 1527 write(msg, '(a,i8,a,a,i8,a,a,a,a,a)' )& 1528 & 'The input value of objbat for atom number ',iatom,ch10,& 1529 & 'is equal to ',objbat(iatom),', an unacceptable value :',ch10,& 1530 & 'it should be between 1 and natom. ',ch10,& 1531 & 'Action: correct objbat in your input file.' 1532 ABI_ERROR(msg) 1533 end if 1534 end do 1535 1536 if(objbn>1)then 1537 do iatom=1,objbn-1 1538 if( objbat(iatom)>=objbat(iatom+1) )then 1539 write(msg, '(a,i8,a,a,a,a,a,a)' )& 1540 & 'The input value of objbat for atom number ',iatom,ch10,& 1541 & 'is larger or equal to the one of the next atom,',ch10,& 1542 & 'while this list should be ordered, and an atom cannot be repeated.',ch10,& 1543 & 'Action: correct the array objbat in the input file.' 1544 ABI_ERROR(msg) 1545 end if 1546 end do 1547 end if 1548 1549 ! Read repetition factors 1550 objbrf(1:3)=1 1551 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objbrf',tread,'INT') 1552 if(tread==1) objbrf(1:3)=intarr(1:3) 1553 write(iout,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:) 1554 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:) 1555 1556 if(tread==1)then 1557 do irep=1,3 1558 if(objbrf(irep)<1)then 1559 write(msg, '(a,a,a,3i8,a,a,a)' )& 1560 & 'The input values of objbrf(1:3) must be positive,',ch10,& 1561 & 'while it is ',objbrf(1:3),'.',ch10,& 1562 & 'Action: correct objbrf in your input file.' 1563 ABI_ERROR(msg) 1564 end if 1565 end do 1566 end if 1567 1568 ! Modify the number of atoms to be read 1569 natom_toberead=natom_toberead-objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1) 1570 ! Read rotations angles and translations 1571 objbro(1:4)=0.0_dp 1572 objbtr(1:12)=0.0_dp 1573 if (objbrf(1)*objbrf(2)*objbrf(3) ==1) then 1574 nread=1 1575 else if (objbrf(2)*objbrf(3) ==1) then 1576 nread=2 1577 else if (objbrf(3) ==1) then 1578 nread=3 1579 else 1580 nread=4 1581 end if 1582 call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objbro',tread,'DPR') 1583 if(tread==1) objbro(1:nread)=dprarr(1:nread) 1584 1585 call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objbtr',tread,'LEN') 1586 if(tread==1) objbtr(1:3*nread)=dprarr(1:3*nread) 1587 1588 write(iout,format01160) 'objbro',objbro(1:4) 1589 write(std_out,format01160) 'objbro',objbro(1:4) 1590 write(iout,format01160) 'objbtr',objbtr(1:12) 1591 write(std_out,format01160) 'objbtr',objbtr(1:12) 1592 1593 ! If needed, read axes, but default to the x-axis to avoid errors later 1594 objbax(1:6)=0.0_dp ; objbax(4)=1.0_dp 1595 if(abs(objbro(1))+abs(objbro(2))+abs(objbro(3))+abs(objbro(4)) > 1.0d-10) then 1596 call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objbax',tread,'LEN') 1597 if(tread==1) objbax(1:6)=dprarr(1:6) 1598 if(tread/=1)then 1599 write(msg, '(a,a,a,a,a,a,a)' )& 1600 & 'The axis of object b (objbax) must be initialized',ch10,& 1601 & 'in the input file, when rotations (objbro) are present.',ch10,& 1602 & 'This is not the case.',ch10,& 1603 & 'Action: initialize objbax in your input file.' 1604 ABI_ERROR(msg) 1605 end if 1606 write(iout,format01160) 'objbax',objbax(1:6) 1607 write(std_out,format01160) 'objbax',objbax(1:6) 1608 end if 1609 axisb(1:3)=objbax(4:6)-objbax(1:3) 1610 normb=axisb(1)**2+axisb(2)**2+axisb(3)**2 1611 if(normb<1.0d-10)then 1612 write(msg, '(5a)' )& 1613 & 'The two points defined by the input array objbax are too',ch10,& 1614 & 'close to each other, and will not be used to define an axis.',ch10,& 1615 & 'Action: correct objbax in your input file.' 1616 ABI_ERROR(msg) 1617 end if 1618 axisb(1:3)=axisb(1:3)/sqrt(normb) 1619 1620 ! Check whether both lists are disjoints. Use a very primitive algorithm. 1621 do iatom=1,objan 1622 do ii=1,objbn 1623 if(objaat(iatom)==objbat(ii))then 1624 write(msg, '(6a,i8,a,i8,3a)' )& 1625 & 'The objects a and b cannot have a common atom, but it is',ch10,& 1626 & 'found that the values of objaat and objbat ',& 1627 & ' are identical, for their',ch10,& 1628 & 'atoms number ',iatom,' and ',ii,'.',ch10,& 1629 & 'Action: change objaat and/or objbat so that they have no common atom anymore.' 1630 ABI_ERROR(msg) 1631 end if 1632 end do 1633 end do 1634 end if ! End condition of existence of a second object 1635 1636 !Check whether the number of atoms to be read obtained by relying 1637 !on natom, vacnum and the object definitions, or from natrd coincide 1638 if(natrd/=natom_toberead)then 1639 write(msg,'(11a,i0,a,i0,2a,i0,a)' )& 1640 & ' ingeobld : ERROR -- ',ch10,& 1641 & ' The number of atoms to be read (natrd) must be equal',ch10,& 1642 & ' to the total number of atoms (natom), plus',ch10,& 1643 & ' the number of vacancies (vacnum), minus',ch10,& 1644 & ' the number of atoms added by the repetition of objects.',ch10,& 1645 & ' This is not the case : natrd= ',natrd,', natom= ',natom,ch10,& 1646 & ', vacnum= ',vacnum,';' 1647 call wrtout(std_out,msg) 1648 1649 if(nobj==1 .or. nobj==2) then 1650 write(msg,'(a,i3,a,3i3,a,i5,a)' )& 1651 & ' object a : objan=',objan,', objarf(1:3)=',objarf(1:3),& 1652 & ' => adds ',objan*(objarf(1)*objarf(2)*objarf(3)-1),' atoms.' 1653 call wrtout(std_out,msg) 1654 end if 1655 1656 if(nobj==2) then 1657 write(msg,'(a,i3,a,3i3,a,i5,a)' )& 1658 & ' object b : objbn=',objbn,', objbrf(1:3)=',objbrf(1:3),& 1659 & ' => adds ',objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1),' atoms.' 1660 call wrtout(std_out,msg) 1661 end if 1662 1663 write(msg,'(3a)' )& 1664 & ' Action : check the correspondence between natom+vacnum on one side,',ch10,& 1665 & ' and natrd, objan, objbn, objarf and objbrf on the other side.' 1666 ABI_ERROR(msg) 1667 end if 1668 1669 !6) Produce full set of atoms 1670 1671 !Print the initial atom coordinates if the geometry builder is used 1672 write(iout, '(/,a)' ) ' Cartesian coordinates of the primitive atoms ' 1673 write(std_out,'(/,a)' )' Cartesian coordinates of the primitive atoms ' 1674 write(iout,format01160) ' ',xcart_read(:,:) 1675 write(std_out,format01160) ' ',xcart_read(:,:) 1676 1677 ABI_MALLOC(typat_full,(natom+vacnum)) 1678 ABI_MALLOC(xcart_full,(3,natom+vacnum)) 1679 1680 !Use the work array xcart_full to produce full set of atoms, 1681 !including those coming from repeated objects. 1682 iatom=1 1683 do iatrd=1,natrd 1684 1685 belonga=0 ; belongb=0 1686 if(nobj==1 .or. nobj==2)then 1687 ! Determine whether the atom belongs to object a 1688 do ii=1,objan 1689 if(iatrd==objaat(ii))belonga=ii 1690 end do 1691 end if 1692 if(nobj==2)then 1693 ! Determine whether the atom belong to object b 1694 do ii=1,objbn 1695 if(iatrd==objbat(ii))belongb=ii 1696 end do 1697 end if 1698 1699 !write(std_out,'(a,i5,a,i2,i2,a)' )' ingeobld : treating iatrd=',iatrd,', belong(a,b)=',belonga,belongb,'.' 1700 1701 ! In case it does not belong to an object 1702 if(belonga==0 .and. belongb==0)then 1703 xcart_full(1:3,iatom)=xcart_read(1:3,iatrd) 1704 typat_full(iatom)=typat_read(iatrd) 1705 iatom=iatom+1 1706 else 1707 1708 ! Repeat, rotate and translate this atom 1709 if(belonga/=0)then 1710 1711 ! Treat object a 1712 ! Compute the relative coordinate of atom with respect to first point of axis 1713 vectora(1:3)=xcart_read(1:3,iatrd)-objaax(1:3) 1714 ! Project on axis 1715 project=vectora(1)*axisa(1)+vectora(2)*axisa(2)+vectora(3)*axisa(3) 1716 ! Get the parallel part 1717 parall(1:3)=project*axisa(1:3) 1718 ! Get the perpendicular part, to be rotated 1719 perpen(1:3)=vectora(1:3)-parall(1:3) 1720 ! Compute the norm of the perpendicular part 1721 norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2 1722 ! Initialisation to avoid warnings even if used behind if rotate == 1. 1723 normper = 0 1724 ! It the norm is too small, there is not need to rotate 1725 rotate=0 1726 if(norm2per>=1.0d-18)then 1727 rotate=1 1728 normper=sqrt(norm2per) 1729 axis2(1:3)=perpen(1:3)/normper 1730 ! Get the vector perpendicular to axisa and axisa2 1731 axis3(1)=axisa(2)*axis2(3)-axisa(3)*axis2(2) 1732 axis3(2)=axisa(3)*axis2(1)-axisa(1)*axis2(3) 1733 axis3(3)=axisa(1)*axis2(2)-axisa(2)*axis2(1) 1734 end if 1735 1736 ! Here the repetition loop 1737 do irep3=1,objarf(3) 1738 do irep2=1,objarf(2) 1739 do irep1=1,objarf(1) 1740 ! Here the rotation 1741 if(rotate==1)then 1742 ! Compute the angle of rotation 1743 angle=objaro(1)+(irep1-1)*objaro(2) + & 1744 & (irep2-1)*objaro(3)+(irep3-1)*objaro(4) 1745 cosine=cos(angle/180.0*pi) 1746 sine=sin(angle/180.0*pi) 1747 rotated(1:3)=objaax(1:3)+parall(1:3)+& 1748 & normper*(cosine*axis2(1:3)+sine*axis3(1:3)) 1749 else 1750 rotated(1:3)=vectora(1:3) 1751 end if 1752 ! Here the translation 1753 xcart_full(1:3,iatom)=rotated(1:3)+objatr(1:3)+& 1754 & (irep1-1)*objatr(4:6)+(irep2-1)*objatr(7:9)+(irep3-1)*objatr(10:12) 1755 typat_full(iatom)=typat_read(iatrd) 1756 iatom=iatom+1 1757 end do 1758 end do 1759 end do ! End the repetition loop 1760 1761 else 1762 ! If the atom belong to object b 1763 ! Compute the relative coordinate of atom with respect to first point of axis 1764 vectorb(1:3)=xcart_read(1:3,iatrd)-objbax(1:3) 1765 ! Project on axis 1766 project=vectorb(1)*axisb(1)+vectorb(2)*axisb(2)+vectorb(3)*axisb(3) 1767 ! Get the parallel part 1768 parall(1:3)=project*axisb(1:3) 1769 ! Get the perpendicular part, to be rotated 1770 perpen(1:3)=vectorb(1:3)-parall(1:3) 1771 ! Compute the norm of the perpendicular part 1772 norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2 1773 ! Initialisation to avoid warnings even if used behind if rotate == 1. 1774 normper = 0 1775 ! It the norm is too small, there is not need to rotate 1776 rotate=0 1777 if(norm2per>=1.0d-18)then 1778 rotate=1 1779 normper=sqrt(norm2per) 1780 axis2(1:3)=perpen(1:3)/normper 1781 ! Get the vector perpendicular to axisb and axis2 1782 axis3(1)=axisb(2)*axis2(3)-axisb(3)*axis2(2) 1783 axis3(2)=axisb(3)*axis2(1)-axisb(1)*axis2(3) 1784 axis3(3)=axisb(1)*axis2(2)-axisb(2)*axis2(1) 1785 end if 1786 ! Here the repetition loop 1787 do irep3=1,objbrf(3) 1788 do irep2=1,objbrf(2) 1789 do irep1=1,objbrf(1) 1790 ! Here the rotation 1791 if(rotate==1)then 1792 ! Compute the angle of rotation 1793 angle=objbro(1)+(irep1-1)*objbro(2) + & 1794 & (irep2-1)*objbro(3)+ (irep3-1)*objbro(4) 1795 cosine=cos(angle/180.0*pi) 1796 sine=sin(angle/180.0*pi) 1797 rotated(1:3)=objbax(1:3)+parall(1:3)+& 1798 & normper*(cosine*axis2(1:3)+sine*axis3(1:3)) 1799 else 1800 rotated(1:3)=vectorb(1:3) 1801 end if 1802 ! Here the translation 1803 xcart_full(1:3,iatom)=rotated(1:3)+objbtr(1:3)+& 1804 & (irep1-1)*objbtr(4:6)+(irep2-1)*objbtr(7:9)+(irep3-1)*objbtr(10:12) 1805 typat_full(iatom)=typat_read(iatrd) 1806 iatom=iatom+1 1807 end do 1808 end do 1809 end do ! End the repetition loop 1810 end if ! Condition of belonging to object b 1811 end if ! Condition of belonging to an object 1812 end do ! Loop on atoms 1813 1814 !Create the vacancies here 1815 if(vacnum/=0)then 1816 ! First label the vacant atoms as belonging to typat 0 1817 do ivac=1,vacnum 1818 typat_full(vaclst(ivac))=0 1819 end do 1820 ! Then compact the arrays 1821 shift=0 1822 do iatom=1,natom 1823 if(typat_full(iatom+shift)==0) shift=shift+1 1824 if(shift/=0)then 1825 xcart_full(1:3,iatom)=xcart_full(1:3,iatom+shift) 1826 typat_full(iatom)=typat_full(iatom+shift) 1827 end if 1828 end do 1829 end if 1830 1831 !Transfer the content of xcart_full and typat_full to the proper location 1832 xcart(:,1:natom)=xcart_full(:,1:natom) 1833 typat(1:natom)=typat_full(1:natom) 1834 1835 ABI_FREE(typat_full) 1836 ABI_FREE(xcart_full) 1837 if(allocated(objaat)) then 1838 ABI_FREE(objaat) 1839 end if 1840 if(allocated(objbat)) then 1841 ABI_FREE(objbat) 1842 end if 1843 1844 ABI_FREE(intarr) 1845 ABI_FREE(dprarr) 1846 if (vacnum>0) then 1847 ABI_FREE(vaclst) 1848 end if 1849 1850 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
2035 subroutine invacuum(jdtset,lenstr,natom,rprimd,string,vacuum,xred) 2036 2037 !Arguments ------------------------------------ 2038 !scalars 2039 integer,intent(in) :: jdtset,lenstr,natom 2040 character(len=*),intent(in) :: string 2041 !arrays 2042 integer,intent(out) :: vacuum(3) 2043 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 2044 2045 !Local variables------------------------------- 2046 !scalars 2047 integer :: ia,ii,marr,tread 2048 real(dp) :: max_diff_xred,ucvol,vacwidth,vacxred 2049 !arrays 2050 integer,allocatable :: list(:) 2051 integer,allocatable :: intarr(:) 2052 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 2053 real(dp),allocatable :: xred_sorted(:) 2054 real(dp),allocatable :: dprarr(:) 2055 2056 ! ************************************************************************* 2057 2058 !Compute the maximum size of arrays intarr and dprarr 2059 marr=3 2060 ABI_MALLOC(intarr,(marr)) 2061 ABI_MALLOC(dprarr,(marr)) 2062 2063 !Get metric quantities 2064 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2065 2066 !Read vacwidth, or set the default 2067 vacwidth=10.0_dp 2068 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacwidth',tread,'LEN') 2069 if(tread==1) vacwidth=dprarr(1) 2070 2071 !Read vacuum, or compute it using the atomic coordinates and vacwidth. 2072 vacuum(1:3)=0 2073 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'vacuum',tread,'INT') 2074 2075 if(tread==1)then 2076 vacuum(1:3)=intarr(1:3) 2077 else 2078 ! For each direction, determine whether a vacuum space exists 2079 ABI_MALLOC(list,(natom)) 2080 ABI_MALLOC(xred_sorted,(natom)) 2081 do ii=1,3 2082 ! This is the minimum xred difference needed to have vacwidth 2083 vacxred=vacwidth*sqrt(sum(gprimd(:,ii)**2)) 2084 ! Project the reduced coordinate in the [0.0_dp,1.0_dp[ interval 2085 xred_sorted(:)=mod(xred(ii,:),1.0_dp) 2086 ! list is dummy 2087 list(:)=0 2088 ! Sort xred_sorted 2089 call sort_dp(natom,xred_sorted,list,tol14) 2090 if(natom==1)then 2091 max_diff_xred=1.0_dp 2092 else 2093 ! Compute the difference between each pair of atom in the sorted order 2094 max_diff_xred=0.0_dp 2095 do ia=1,natom-1 2096 max_diff_xred=max(max_diff_xred,xred_sorted(ia+1)-xred_sorted(ia)) 2097 end do 2098 ! Do not forget the image of the first atom in the next cell 2099 max_diff_xred=max(max_diff_xred,1.0_dp+xred_sorted(1)-xred_sorted(ia)) 2100 end if 2101 if(vacxred<max_diff_xred+tol10)vacuum(ii)=1 2102 end do 2103 ABI_FREE(list) 2104 ABI_FREE(xred_sorted) 2105 end if 2106 2107 !DEBUG 2108 !write(std_out,*)' invacuum : vacuum=',vacuum(1:3) 2109 !ENDDEBUG 2110 2111 ABI_FREE(intarr) 2112 ABI_FREE(dprarr) 2113 2114 end subroutine invacuum