TABLE OF CONTENTS
ABINIT/ingeobld [ 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.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG) 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
iout=unit number of output file jdtset=number of the dataset looked for lenstr=actual length of the string natrd=number of atoms that have been read in the calling routine natom=number of atoms nobj=the number of objects string*(*)=character string containing all the input data, used only if choice=1 or 3. Initialized previously in instrng. typat_read(natrd)=type integer for each atom in the primitive set xcart_read(3,natrd)=cartesian coordinates of atoms (bohr), in the primitive set
OUTPUT
typat(natom)=type integer for each atom in cell xcart(3,natom)=cartesian coordinates of atoms (bohr)
PARENTS
ingeo
CHILDREN
intagm,wrtout
SOURCE
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 #include "abi_common.h" 49 50 subroutine ingeobld (iout,jdtset,lenstr,natrd,natom,nobj,string,typat,typat_read,xcart,xcart_read) 51 52 use defs_basis 53 use m_profiling_abi 54 use m_errors 55 56 use m_parser, only : intagm 57 58 !This section has been created automatically by the script Abilint (TD). 59 !Do not modify the following lines by hand. 60 #undef ABI_FUNC 61 #define ABI_FUNC 'ingeobld' 62 use interfaces_14_hidewrite 63 !End of the abilint section 64 65 implicit none 66 67 !Arguments ------------------------------------ 68 !scalars 69 integer,intent(in) :: iout,jdtset,lenstr,natom,natrd,nobj 70 character(len=*),intent(in) :: string 71 !arrays 72 integer,intent(in) :: typat_read(natrd) 73 integer,intent(out) :: typat(natom) 74 real(dp),intent(in) :: xcart_read(3,natrd) 75 real(dp),intent(out) :: xcart(3,natom) 76 77 !Local variables------------------------------- 78 character(len=*), parameter :: format01110 ="(1x,a6,1x,(t9,8i8) )" 79 character(len=*), parameter :: format01160 ="(1x,a6,1x,1p,(t9,3g18.10)) " 80 !scalars 81 integer :: belonga,belongb,iatom,iatrd,ii,irep,irep1,irep2,irep3,ivac,marr 82 integer :: natom_toberead,nread,objan,objbn,rotate,shift,tread,vacnum 83 real(dp) :: angle,cosine,norm2per,norma,normb,normper,project,sine 84 character(len=500) :: message 85 !arrays 86 integer :: objarf(3),objbrf(3) 87 integer,allocatable :: objaat(:),objbat(:),vaclst(:) 88 real(dp) :: axis2(3),axis3(3),axisa(3),axisb(3),objaax(6),objaro(4),objatr(12) 89 real(dp) :: objbax(6),objbro(4),objbtr(12),parall(3),perpen(3),rotated(3) 90 real(dp) :: vectora(3),vectorb(3) 91 real(dp),allocatable :: typat_full(:),xcart_full(:,:) 92 !no_abirules 93 !Dummy arguments for subroutine 'intagm' to parse input file 94 integer,allocatable :: intarr(:) 95 real(dp),allocatable :: dprarr(:) 96 97 ! ************************************************************************* 98 99 marr=max(12,3*natom) 100 ABI_ALLOCATE(intarr,(marr)) 101 ABI_ALLOCATE(dprarr,(marr)) 102 103 !1) Set up the number of vacancies. 104 105 !This is the default 106 vacnum=0 107 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacnum',tread,'INT') 108 if(tread==1) vacnum=intarr(1) 109 110 if (vacnum>0)then 111 ABI_ALLOCATE(vaclst,(vacnum)) 112 ! Read list of atoms to be suppressed to create vacancies 113 call intagm(dprarr,intarr,jdtset,marr,vacnum,string(1:lenstr),'vaclst',tread,'INT') 114 if(tread==1) vaclst(:)=intarr(1:vacnum) 115 if(tread/=1)then 116 write(message, '(a,a,a,a,a)' )& 117 & 'The array vaclst MUST be initialized in the input file',ch10,& 118 & 'when vacnum is non-zero.',ch10,& 119 & 'Action: initialize vaclst in your input file.' 120 MSG_ERROR(message) 121 end if 122 end if 123 124 natom_toberead=natom+vacnum 125 126 !2) Set up list and number of atoms in objects, and the -------------- 127 !operations to be performed on objects. 128 129 write(message,'(80a,a)')('=',ii=1,80),ch10 130 call wrtout(std_out,message,'COLL') 131 call wrtout(iout,message,'COLL') 132 133 write(message, '(a,a)' )& 134 & '--ingeobld: echo values of variables connected to objects --------',ch10 135 call wrtout(std_out,message,'COLL') 136 call wrtout(iout,message,'COLL') 137 138 if(vacnum>0)then 139 write(iout,format01110) 'vacnum',vacnum 140 write(std_out,format01110) 'vacnum',vacnum 141 write(iout,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:) 142 write(std_out,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:) 143 write(iout, '(a)' ) ' ' 144 write(std_out,'(a)' ) ' ' 145 end if 146 147 write(iout,format01110) 'nobj',nobj 148 write(std_out,format01110) 'nobj',nobj 149 150 if(nobj/=1 .and. nobj/=2)then 151 write(message, '(a,a,a,i8,a,a,a)' )& 152 & 'The number of object (nobj) must be either 1 or 2,',ch10,& 153 & 'while the input file has nobj=',nobj,'.',ch10,& 154 & 'Action: correct nobj in your input file.' 155 MSG_ERROR(message) 156 end if 157 158 if(nobj==1 .or. nobj==2)then 159 160 ! Read the number of atoms of the object a 161 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objan',tread,'INT') 162 if(tread==1) objan=intarr(1) 163 164 if(tread/=1)then 165 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 166 & 'The number of atoms in object a (objan) must be initialized',ch10,& 167 & 'in the input file, when nobj=',nobj,'.',ch10,& 168 & 'This is not the case.',ch10,& 169 & 'Action: correct objan in your input file.' 170 MSG_ERROR(message) 171 end if 172 173 write(iout, '(a)' ) ' ' 174 write(std_out,'(a)' ) ' ' 175 write(iout,format01110) 'objan',objan 176 write(std_out,format01110) 'objan',objan 177 178 if(objan<=1 .or. objan>natom)then 179 write(message, '(a,a,a,a,a,i8,a,a,a)' )& 180 & 'The number of atoms in object a (objan) must be larger than 0',ch10,& 181 & 'and smaller than natom.',ch10,& 182 & 'It is equal to ',objan,', an unacceptable value.',ch10,& 183 & 'Action: correct objan in your input file.' 184 MSG_ERROR(message) 185 end if 186 187 ! Read list of atoms in object a 188 call intagm(dprarr,intarr,jdtset,marr,objan,string(1:lenstr),'objaat',tread,'INT') 189 ABI_ALLOCATE(objaat,(objan)) 190 if(tread==1) objaat(1:objan)=intarr(1:objan) 191 192 if(tread/=1)then 193 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 194 & 'The list of atoms in object a (objaat) must be initialized',ch10,& 195 & 'in the input file, when nobj=',nobj,'.',ch10,& 196 & 'This is not the case.',ch10,& 197 & 'Action: initialize objaat in your input file.' 198 MSG_ERROR(message) 199 end if 200 201 write(iout,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:) 202 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:) 203 204 do iatom=1,objan 205 if(objaat(iatom)<1 .or. objaat(iatom)>natom)then 206 write(message, '(a,i8,a,a,i8,4a)' )& 207 & 'The input value of objaat for atom number ',iatom,ch10,& 208 & 'is equal to ',objaat(iatom),', an unacceptable value :',ch10,& 209 & 'it should be between 1 and natom. ',& 210 & 'Action: correct the array objaat in your input file.' 211 MSG_ERROR(message) 212 end if 213 end do 214 215 if(objan>1)then 216 do iatom=1,objan-1 217 if( objaat(iatom)>=objaat(iatom+1) )then 218 write(message, '(a,i8,a,a,a,a,a,a)' )& 219 & 'The input value of objaat for atom number ',iatom,ch10,& 220 & 'is larger or equal to the one of the next atom,',ch10,& 221 & 'while this list should be ordered, and an atom cannot be repeated.',ch10,& 222 & 'Action: correct the array objaat in your input file.' 223 MSG_ERROR(message) 224 end if 225 end do 226 end if 227 228 ! Read repetition factors 229 objarf(1:3)=1 230 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objarf',tread,'INT') 231 if(tread==1) objarf(1:3)=intarr(1:3) 232 write(iout,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:) 233 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:) 234 235 if(tread==1)then 236 do irep=1,3 237 if(objarf(irep)<1)then 238 write(message, '(a,a,a,3i8,a,a,a)' )& 239 & 'The input values of objarf(1:3) must be positive,',ch10,& 240 & 'while it is ',objarf(1:3),'.',ch10,& 241 & 'Action: correct objarf in your input file.' 242 MSG_ERROR(message) 243 end if 244 end do 245 end if 246 247 ! Modify the number of atoms to be read 248 natom_toberead=natom_toberead-objan*(objarf(1)*objarf(2)*objarf(3)-1) 249 250 ! Read rotations angles and translations 251 objaro(1:4)=0.0_dp 252 objatr(1:12)=0.0_dp 253 if (objarf(1)*objarf(2)*objarf(3) ==1) then 254 nread=1 255 else if (objarf(2)*objarf(3) ==1) then 256 nread=2 257 else if (objarf(3) ==1) then 258 nread=3 259 else 260 nread=4 261 end if 262 call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objaro',tread,'DPR') 263 if(tread==1) objaro(1:nread)=dprarr(1:nread) 264 265 call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objatr',tread,'LEN') 266 267 if(tread==1) objatr(1:3*nread)=dprarr(1:3*nread) 268 write(iout,format01160) 'objaro',objaro(1:4) 269 write(std_out,format01160) 'objaro',objaro(1:4) 270 write(iout,format01160) 'objatr',objatr(1:12) 271 write(std_out,format01160) 'objatr',objatr(1:12) 272 ! If needed, read axes, but default to the x-axis to avoid errors later 273 objaax(1:6)=0.0_dp ; objaax(4)=1.0_dp 274 275 if(abs(objaro(1))+abs(objaro(2))+abs(objaro(3))+abs(objaro(4)) > 1.0d-10) then 276 call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objaax',tread,'LEN') 277 if(tread==1) objaax(1:6)=dprarr(1:6) 278 if(tread/=1)then 279 write(message, '(a,a,a,a,a,a,a)' )& 280 & 'The axis of object a (objaax) must be initialized',ch10,& 281 & 'in the input file, when rotations (objaro) are present.',ch10,& 282 & 'This is not the case.',ch10,& 283 & 'Action: initialize objaax in your input file.' 284 MSG_ERROR(message) 285 end if 286 write(iout,format01160) 'objaax',objaax(1:6) 287 write(std_out,format01160) 'objaax',objaax(1:6) 288 end if 289 290 axisa(1:3)=objaax(4:6)-objaax(1:3) 291 norma=axisa(1)**2+axisa(2)**2+axisa(3)**2 292 293 if(norma<1.0d-10)then 294 write(message, '(5a)' )& 295 & 'The two points defined by the input array objaax are too',ch10,& 296 & 'close to each other, and will not be used to define an axis.',ch10,& 297 & 'Action: correct objaax in your input file.' 298 MSG_ERROR(message) 299 end if 300 axisa(1:3)=axisa(1:3)/sqrt(norma) 301 302 ! End condition of existence of a first object 303 end if 304 305 if(nobj==2)then 306 307 ! Read the number of atoms of the object b 308 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objbn',tread,'INT') 309 if(tread==1) objbn=intarr(1) 310 311 if(tread/=1)then 312 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 313 & 'The number of atoms in object b (objbn) must be initialized',ch10,& 314 & 'in the input file, when nobj=',nobj,'.',ch10,& 315 & 'This is not the case.',ch10,& 316 & 'Action: initialize objbn in your input file.' 317 MSG_ERROR(message) 318 end if 319 320 write(iout, '(a)' ) ' ' 321 write(std_out,'(a)' ) ' ' 322 write(iout,format01110) 'objbn',objbn 323 write(std_out,format01110) 'objbn',objbn 324 325 if(objbn<=1 .or. objbn>natom)then 326 write(message, '(a,a,a,a,a,i8,a,a,a)' )& 327 & 'The number of atoms in object b (objbn) must be larger than 0',ch10,& 328 & 'and smaller than natom.',ch10,& 329 & 'It is equal to ',objbn,', an unacceptable value.',ch10,& 330 & 'Action: correct objbn in your input file.' 331 MSG_ERROR(message) 332 end if 333 334 ! Read list of atoms in object b 335 call intagm(dprarr,intarr,jdtset,marr,objbn,string(1:lenstr),'objbat',tread,'INT') 336 ABI_ALLOCATE(objbat,(objbn)) 337 338 if(tread==1) objbat(1:objbn)=intarr(1:objbn) 339 if(tread/=1)then 340 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 341 & 'The list of atoms in object b (objbat) must be initialized',ch10,& 342 & 'in the input file, when nobj=',nobj,'.',ch10,& 343 & 'This is not the case.',ch10,& 344 & 'Action: initialize objbat in your input file.' 345 MSG_ERROR(message) 346 end if 347 348 write(iout,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:) 349 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:) 350 351 do iatom=1,objbn 352 if(objbat(iatom)<1 .or. objbat(iatom)>natom)then 353 write(message, '(a,i8,a,a,i8,a,a,a,a,a)' )& 354 & 'The input value of objbat for atom number ',iatom,ch10,& 355 & 'is equal to ',objbat(iatom),', an unacceptable value :',ch10,& 356 & 'it should be between 1 and natom. ',ch10,& 357 & 'Action: correct objbat in your input file.' 358 MSG_ERROR(message) 359 end if 360 end do 361 362 if(objbn>1)then 363 do iatom=1,objbn-1 364 if( objbat(iatom)>=objbat(iatom+1) )then 365 write(message, '(a,i8,a,a,a,a,a,a)' )& 366 & 'The input value of objbat for atom number ',iatom,ch10,& 367 & 'is larger or equal to the one of the next atom,',ch10,& 368 & 'while this list should be ordered, and an atom cannot be repeated.',ch10,& 369 & 'Action: correct the array objbat in the input file.' 370 MSG_ERROR(message) 371 end if 372 end do 373 end if 374 375 ! Read repetition factors 376 objbrf(1:3)=1 377 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objbrf',tread,'INT') 378 if(tread==1) objbrf(1:3)=intarr(1:3) 379 write(iout,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:) 380 write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:) 381 382 if(tread==1)then 383 do irep=1,3 384 if(objbrf(irep)<1)then 385 write(message, '(a,a,a,3i8,a,a,a)' )& 386 & 'The input values of objbrf(1:3) must be positive,',ch10,& 387 & 'while it is ',objbrf(1:3),'.',ch10,& 388 & 'Action: correct objbrf in your input file.' 389 MSG_ERROR(message) 390 end if 391 end do 392 end if 393 394 ! Modify the number of atoms to be read 395 natom_toberead=natom_toberead-objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1) 396 ! Read rotations angles and translations 397 objbro(1:4)=0.0_dp 398 objbtr(1:12)=0.0_dp 399 if (objbrf(1)*objbrf(2)*objbrf(3) ==1) then 400 nread=1 401 else if (objbrf(2)*objbrf(3) ==1) then 402 nread=2 403 else if (objbrf(3) ==1) then 404 nread=3 405 else 406 nread=4 407 end if 408 call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objbro',tread,'DPR') 409 if(tread==1) objbro(1:nread)=dprarr(1:nread) 410 411 call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objbtr',tread,'LEN') 412 if(tread==1) objbtr(1:3*nread)=dprarr(1:3*nread) 413 414 write(iout,format01160) 'objbro',objbro(1:4) 415 write(std_out,format01160) 'objbro',objbro(1:4) 416 write(iout,format01160) 'objbtr',objbtr(1:12) 417 write(std_out,format01160) 'objbtr',objbtr(1:12) 418 419 ! If needed, read axes, but default to the x-axis to avoid errors later 420 objbax(1:6)=0.0_dp ; objbax(4)=1.0_dp 421 if(abs(objbro(1))+abs(objbro(2))+abs(objbro(3))+abs(objbro(4)) > 1.0d-10) then 422 call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objbax',tread,'LEN') 423 if(tread==1) objbax(1:6)=dprarr(1:6) 424 if(tread/=1)then 425 write(message, '(a,a,a,a,a,a,a)' )& 426 & 'The axis of object b (objbax) must be initialized',ch10,& 427 & 'in the input file, when rotations (objbro) are present.',ch10,& 428 & 'This is not the case.',ch10,& 429 & 'Action: initialize objbax in your input file.' 430 MSG_ERROR(message) 431 end if 432 write(iout,format01160) 'objbax',objbax(1:6) 433 write(std_out,format01160) 'objbax',objbax(1:6) 434 end if 435 axisb(1:3)=objbax(4:6)-objbax(1:3) 436 normb=axisb(1)**2+axisb(2)**2+axisb(3)**2 437 if(normb<1.0d-10)then 438 write(message, '(5a)' )& 439 & 'The two points defined by the input array objbax are too',ch10,& 440 & 'close to each other, and will not be used to define an axis.',ch10,& 441 & 'Action: correct objbax in your input file.' 442 MSG_ERROR(message) 443 end if 444 axisb(1:3)=axisb(1:3)/sqrt(normb) 445 446 ! Check whether both lists are disjoints. Use a very primitive algorithm. 447 do iatom=1,objan 448 do ii=1,objbn 449 if(objaat(iatom)==objbat(ii))then 450 write(message, '(6a,i8,a,i8,3a)' )& 451 & 'The objects a and b cannot have a common atom, but it is',ch10,& 452 & 'found that the values of objaat and objbat ',& 453 & ' are identical, for their',ch10,& 454 & 'atoms number ',iatom,' and ',ii,'.',ch10,& 455 & 'Action: change objaat and/or objbat so that they have no common atom anymore.' 456 MSG_ERROR(message) 457 end if 458 end do 459 end do 460 461 ! End condition of existence of a second object 462 end if 463 464 !Check whether the number of atoms to be read obtained by relying 465 !on natom, vacnum and the object definitions, or from natrd coincide 466 if(natrd/=natom_toberead)then 467 write(message,'(11a,i0,a,i0,2a,i0,a)' )& 468 & ' ingeobld : ERROR -- ',ch10,& 469 & ' The number of atoms to be read (natrd) must be equal',ch10,& 470 & ' to the total number of atoms (natom), plus',ch10,& 471 & ' the number of vacancies (vacnum), minus',ch10,& 472 & ' the number of atoms added by the repetition of objects.',ch10,& 473 & ' This is not the case : natrd= ',natrd,', natom= ',natom,ch10,& 474 & ', vacnum= ',vacnum,';' 475 call wrtout(std_out,message,"COLL") 476 477 if(nobj==1 .or. nobj==2) then 478 write(message,'(a,i3,a,3i3,a,i5,a)' )& 479 & ' object a : objan=',objan,', objarf(1:3)=',objarf(1:3),& 480 & ' => adds ',objan*(objarf(1)*objarf(2)*objarf(3)-1),' atoms.' 481 call wrtout(std_out,message,"COLL") 482 end if 483 484 if(nobj==2) then 485 write(message,'(a,i3,a,3i3,a,i5,a)' )& 486 & ' object b : objbn=',objbn,', objbrf(1:3)=',objbrf(1:3),& 487 & ' => adds ',objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1),' atoms.' 488 call wrtout(std_out,message,"COLL") 489 end if 490 491 write(message,'(3a)' )& 492 & ' Action : check the correspondence between natom+vacnum on one side,',ch10,& 493 & ' and natrd, objan, objbn, objarf and objbrf on the other side.' 494 MSG_ERROR(message) 495 end if 496 497 !6) Produce full set of atoms 498 499 !Print the initial atom coordinates if the geometry builder is used 500 write(iout, '(/,a)' ) ' Cartesian coordinates of the primitive atoms ' 501 write(std_out,'(/,a)' )' Cartesian coordinates of the primitive atoms ' 502 write(iout,format01160) ' ',xcart_read(:,:) 503 write(std_out,format01160) ' ',xcart_read(:,:) 504 505 ABI_ALLOCATE(typat_full,(natom+vacnum)) 506 ABI_ALLOCATE(xcart_full,(3,natom+vacnum)) 507 508 !Use the work array xcart_full to produce full set of atoms, 509 !including those coming from repeated objects. 510 iatom=1 511 do iatrd=1,natrd 512 513 belonga=0 ; belongb=0 514 if(nobj==1 .or. nobj==2)then 515 ! Determine whether the atom belongs to object a 516 do ii=1,objan 517 if(iatrd==objaat(ii))belonga=ii 518 end do 519 end if 520 if(nobj==2)then 521 ! Determine whether the atom belong to object b 522 do ii=1,objbn 523 if(iatrd==objbat(ii))belongb=ii 524 end do 525 end if 526 527 write(std_out,'(a,i5,a,i2,i2,a)' ) & 528 & ' ingeobld : treating iatrd=',iatrd,', belong(a,b)=',belonga,belongb,'.' 529 530 ! In case it does not belong to an object 531 if(belonga==0 .and. belongb==0)then 532 xcart_full(1:3,iatom)=xcart_read(1:3,iatrd) 533 typat_full(iatom)=typat_read(iatrd) 534 iatom=iatom+1 535 else 536 537 ! Repeat, rotate and translate this atom 538 if(belonga/=0)then 539 540 ! Treat object a 541 ! Compute the relative coordinate of atom with respect to first point of axis 542 vectora(1:3)=xcart_read(1:3,iatrd)-objaax(1:3) 543 ! Project on axis 544 project=vectora(1)*axisa(1)+vectora(2)*axisa(2)+vectora(3)*axisa(3) 545 ! Get the parallel part 546 parall(1:3)=project*axisa(1:3) 547 ! Get the perpendicular part, to be rotated 548 perpen(1:3)=vectora(1:3)-parall(1:3) 549 ! Compute the norm of the perpendicular part 550 norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2 551 ! Initialisation to avoid warnings even if used behind if rotate == 1. 552 normper = 0 553 ! It the norm is too small, there is not need to rotate 554 rotate=0 555 if(norm2per>=1.0d-18)then 556 rotate=1 557 normper=sqrt(norm2per) 558 axis2(1:3)=perpen(1:3)/normper 559 ! Get the vector perpendicular to axisa and axisa2 560 axis3(1)=axisa(2)*axis2(3)-axisa(3)*axis2(2) 561 axis3(2)=axisa(3)*axis2(1)-axisa(1)*axis2(3) 562 axis3(3)=axisa(1)*axis2(2)-axisa(2)*axis2(1) 563 end if 564 565 ! Here the repetition loop 566 do irep3=1,objarf(3) 567 do irep2=1,objarf(2) 568 do irep1=1,objarf(1) 569 ! Here the rotation 570 if(rotate==1)then 571 ! Compute the angle of rotation 572 angle=objaro(1)+(irep1-1)*objaro(2)+ & 573 & (irep2-1)*objaro(3)+(irep3-1)*objaro(4) 574 cosine=cos(angle/180.0*pi) 575 sine=sin(angle/180.0*pi) 576 rotated(1:3)=objaax(1:3)+parall(1:3)+& 577 & normper*(cosine*axis2(1:3)+sine*axis3(1:3)) 578 else 579 rotated(1:3)=vectora(1:3) 580 end if 581 ! Here the translation 582 xcart_full(1:3,iatom)=rotated(1:3)+objatr(1:3)+& 583 & (irep1-1)*objatr(4:6)+(irep2-1)*objatr(7:9)+(irep3-1)*objatr(10:12) 584 typat_full(iatom)=typat_read(iatrd) 585 iatom=iatom+1 586 end do 587 end do 588 ! End the repetition loop 589 end do 590 591 else 592 ! If the atom belong to object b 593 ! Compute the relative coordinate of atom with respect to first point of axis 594 vectorb(1:3)=xcart_read(1:3,iatrd)-objbax(1:3) 595 ! Project on axis 596 project=vectorb(1)*axisb(1)+vectorb(2)*axisb(2)+vectorb(3)*axisb(3) 597 ! Get the parallel part 598 parall(1:3)=project*axisb(1:3) 599 ! Get the perpendicular part, to be rotated 600 perpen(1:3)=vectorb(1:3)-parall(1:3) 601 ! Compute the norm of the perpendicular part 602 norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2 603 ! Initialisation to avoid warnings even if used behind if rotate == 1. 604 normper = 0 605 ! It the norm is too small, there is not need to rotate 606 rotate=0 607 if(norm2per>=1.0d-18)then 608 rotate=1 609 normper=sqrt(norm2per) 610 axis2(1:3)=perpen(1:3)/normper 611 ! Get the vector perpendicular to axisb and axis2 612 axis3(1)=axisb(2)*axis2(3)-axisb(3)*axis2(2) 613 axis3(2)=axisb(3)*axis2(1)-axisb(1)*axis2(3) 614 axis3(3)=axisb(1)*axis2(2)-axisb(2)*axis2(1) 615 end if 616 ! Here the repetition loop 617 do irep3=1,objbrf(3) 618 do irep2=1,objbrf(2) 619 do irep1=1,objbrf(1) 620 ! Here the rotation 621 if(rotate==1)then 622 ! Compute the angle of rotation 623 angle=objbro(1)+(irep1-1)*objbro(2)+ & 624 & (irep2-1)*objbro(3)+ (irep3-1)*objbro(4) 625 cosine=cos(angle/180.0*pi) 626 sine=sin(angle/180.0*pi) 627 rotated(1:3)=objbax(1:3)+parall(1:3)+& 628 & normper*(cosine*axis2(1:3)+sine*axis3(1:3)) 629 else 630 rotated(1:3)=vectorb(1:3) 631 end if 632 ! Here the translation 633 xcart_full(1:3,iatom)=rotated(1:3)+objbtr(1:3)+& 634 & (irep1-1)*objbtr(4:6)+(irep2-1)*objbtr(7:9)+(irep3-1)*objbtr(10:12) 635 typat_full(iatom)=typat_read(iatrd) 636 iatom=iatom+1 637 end do 638 end do 639 ! End the repetition loop 640 end do 641 642 ! End the condition of belonging to object b 643 end if 644 645 ! End the condition of belonging to an object 646 end if 647 648 ! End the loop on atoms 649 end do 650 651 !Create the vacancies here 652 if(vacnum/=0)then 653 ! First label the vacant atoms as belonging to typat 0 654 do ivac=1,vacnum 655 typat_full(vaclst(ivac))=0 656 end do 657 ! Then compact the arrays 658 shift=0 659 do iatom=1,natom 660 if(typat_full(iatom+shift)==0) shift=shift+1 661 if(shift/=0)then 662 xcart_full(1:3,iatom)=xcart_full(1:3,iatom+shift) 663 typat_full(iatom)=typat_full(iatom+shift) 664 end if 665 end do 666 end if 667 668 !Transfer the content of xcart_full and typat_full to the proper 669 !location 670 xcart(:,1:natom)=xcart_full(:,1:natom) 671 typat(1:natom)=typat_full(1:natom) 672 673 ABI_DEALLOCATE(typat_full) 674 ABI_DEALLOCATE(xcart_full) 675 if(allocated(objaat)) then 676 ABI_DEALLOCATE(objaat) 677 end if 678 if(allocated(objbat)) then 679 ABI_DEALLOCATE(objbat) 680 end if 681 682 ABI_DEALLOCATE(intarr) 683 ABI_DEALLOCATE(dprarr) 684 if (vacnum>0) then 685 ABI_DEALLOCATE(vaclst) 686 end if 687 688 end subroutine ingeobld