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