1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch 2!! Copyright (C) 2013 M. Oliveira 3!! 4!! This program is free software; you can redistribute it and/or modify 5!! it under the terms of the GNU General Public License as published by 6!! the Free Software Foundation; either version 2, or (at your option) 7!! any later version. 8!! 9!! This program is distributed in the hope that it will be useful, 10!! but WITHOUT ANY WARRANTY; without even the implied warranty of 11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12!! GNU General Public License for more details. 13!! 14!! You should have received a copy of the GNU General Public License 15!! along with this program; if not, write to the Free Software 16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 17!! 02110-1301, USA. 18!! 19 20#include "global.h" 21 22module simul_box_oct_m 23 use atom_oct_m 24 use iso_c_binding 25 use gdlib_oct_m 26 use geometry_oct_m 27 use global_oct_m 28 use io_oct_m 29 use kpoints_oct_m 30 use lalg_basic_oct_m 31 use lookup_oct_m 32 use math_oct_m 33 use messages_oct_m 34 use mpi_oct_m 35 use namespace_oct_m 36 use parser_oct_m 37 use profiling_oct_m 38 use space_oct_m 39 use species_oct_m 40 use string_oct_m 41 use symm_op_oct_m 42 use symmetries_oct_m 43 use unit_oct_m 44 use unit_system_oct_m 45 use varinfo_oct_m 46 47 implicit none 48 49 private 50 public :: & 51 simul_box_t, & 52 simul_box_init, & 53 simul_box_lookup_init, & 54 simul_box_interp_init, & 55 simul_box_dump, & 56 simul_box_load, & 57 simul_box_end, & 58 simul_box_write_info, & 59 simul_box_write_short_info, & 60 simul_box_is_periodic, & 61 simul_box_has_zero_bc, & 62 simul_box_in_box, & 63 simul_box_in_box_vec, & 64 simul_box_atoms_in_box, & 65 simul_box_copy, & 66 simul_box_periodic_atom_in_box, & 67 simul_box_symmetry_check, & 68 reciprocal_lattice, & 69 interp_t, & 70 multiresolution_t 71 72 integer, parameter, public :: & 73 SPHERE = 1, & 74 CYLINDER = 2, & 75 MINIMUM = 3, & 76 PARALLELEPIPED = 4, & 77 BOX_IMAGE = 5, & 78 HYPERCUBE = 6, & 79 BOX_USDEF = 77 80 !< BOX_USDEF shares a number with other 'user_defined' input file options. 81 82 type :: interp_t 83 ! Components are public by default 84 integer :: nn, order !< interpolation points and order 85 FLOAT, pointer :: ww(:) !< weights 86 integer, pointer :: posi(:) !< positions 87 end type interp_t 88 89 90 type :: multiresolution_t 91 ! Components are public by default 92 type(interp_t) :: interp !< interpolation points 93 integer, private :: num_areas !< number of multiresolution areas 94 integer :: num_radii !< number of radii (resolution borders) 95 FLOAT, pointer :: radius(:) !< radius of the high-resolution area 96 FLOAT :: center(MAX_DIM) !< central point 97 end type multiresolution_t 98 99 type simul_box_t 100 ! Components are public by default 101 type(symmetries_t) :: symm 102 !> 1->sphere, 2->cylinder, 3->sphere around each atom, 103 !! 4->parallelepiped (orthonormal, up to now). 104 integer :: box_shape 105 106 FLOAT :: rsize !< the radius of the sphere or of the cylinder 107 FLOAT :: xsize !< the length of the cylinder in the x-direction 108 FLOAT :: lsize(MAX_DIM) !< half of the length of the parallelepiped in each direction. 109 110 type(lookup_t), private :: atom_lookup 111 112 character(len=1024), private :: user_def !< for the user-defined box 113 114 logical :: mr_flag !< .true. when using multiresolution 115 type(multiresolution_t) :: hr_area !< high-resolution areas 116 117 FLOAT :: rlattice_primitive(MAX_DIM,MAX_DIM) !< lattice primitive vectors 118 FLOAT :: rlattice (MAX_DIM,MAX_DIM) !< lattice vectors 119 FLOAT :: klattice_primitive(MAX_DIM,MAX_DIM) !< reciprocal-lattice primitive vectors 120 FLOAT :: klattice (MAX_DIM,MAX_DIM) !< reciprocal-lattice vectors 121 FLOAT, private :: volume_element !< the volume element in real space 122 FLOAT :: surface_element (MAX_DIM) !< surface element in real space 123 FLOAT :: rcell_volume !< the volume of the cell in real space 124 FLOAT :: alpha, beta, gamma !< the angles defining the cell 125 FLOAT :: rmetric (MAX_DIM,MAX_DIM) !< metric for the real space lattice vectors 126 FLOAT :: stress_tensor(MAX_DIM,MAX_DIM) !< reciprocal-lattice primitive vectors 127 logical :: nonorthogonal 128 129 type(kpoints_t) :: kpoints !< the k-points 130 131 integer :: dim 132 integer :: periodic_dim 133 134 !> for the box defined through an image 135 integer :: image_size(1:2) 136 type(c_ptr), private :: image 137 character(len=200), private :: filename 138 139 end type simul_box_t 140 141 character(len=22), parameter :: dump_tag = '*** simul_box_dump ***' 142 143contains 144 145 !-------------------------------------------------------------- 146 subroutine simul_box_init(sb, namespace, geo, space) 147 type(simul_box_t), intent(inout) :: sb 148 type(namespace_t), intent(in) :: namespace 149 type(geometry_t), intent(inout) :: geo 150 type(space_t), intent(in) :: space 151 152 ! some local stuff 153 FLOAT :: def_h, def_rsize 154 logical :: only_gamma_kpoint 155 156 PUSH_SUB(simul_box_init) 157 158 call geometry_grid_defaults(geo, def_h, def_rsize) 159 160 call read_misc() ! Miscellaneous stuff. 161 call read_box() ! Parameters defining the simulation box. 162 call simul_box_lookup_init(sb, geo) 163 call simul_box_build_lattice(sb, namespace) ! Build lattice vectors. 164 call simul_box_atoms_in_box(sb, geo, namespace, .true.) ! Put all the atoms inside the box. 165 166 call simul_box_check_atoms_are_too_close(geo, sb, namespace) 167 168 call symmetries_init(sb%symm, namespace, geo, sb%dim, sb%periodic_dim, sb%rlattice, sb%klattice) 169 170 ! we need k-points for periodic systems 171 only_gamma_kpoint = (sb%periodic_dim == 0) 172 call kpoints_init(sb%kpoints, namespace, sb%symm, sb%dim, sb%rlattice, sb%klattice, only_gamma_kpoint) 173 174 call simul_box_symmetry_check(sb, geo, sb%dim, namespace) 175 176 POP_SUB(simul_box_init) 177 178 contains 179 180 181 !-------------------------------------------------------------- 182 subroutine read_misc() 183 184 integer :: idir, irad, order 185 type(block_t) :: blk 186 187 PUSH_SUB(simul_box_init.read_misc) 188 189 sb%dim = space%dim 190 191 !%Variable PeriodicDimensions 192 !%Type integer 193 !%Default 0 194 !%Section System 195 !%Description 196 !% Define how many directions are to be considered periodic. It has to be a number 197 !% between zero and <tt>Dimensions</tt>. 198 !% (WARNING: For systems that are periodic in 1D and 2D, interaction between ions is assumed to be periodic in 3D. 199 !% This affects the calculation of total energy and forces.) 200 !%Option 0 201 !% No direction is periodic (molecule). 202 !%Option 1 203 !% The <i>x</i> direction is periodic (wire, polymer). 204 !%Option 2 205 !% The <i>x</i> and <i>y</i> directions are periodic (slab). 206 !%Option 3 207 !% The <i>x</i>, <i>y</i>, and <i>z</i> directions are periodic (bulk). 208 !%End 209 210 if(geo%periodic_dim == -1) then 211 call parse_variable(namespace, 'PeriodicDimensions', 0, sb%periodic_dim) 212 else 213 sb%periodic_dim = geo%periodic_dim 214 end if 215 if ((sb%periodic_dim < 0) .or. (sb%periodic_dim > MAX_DIM) .or. (sb%periodic_dim > sb%dim)) then 216 call messages_input_error(namespace, 'PeriodicDimensions') 217 end if 218 219 if(sb%periodic_dim > 0 .and. sb%periodic_dim < sb%dim) then 220 call messages_experimental('Support for mixed periodicity systems') 221 end if 222 223 if(sb%periodic_dim == 1) then 224 call messages_write('For systems that are periodic in 1D, interaction between', new_line = .true.) 225 call messages_write('ions is assumed to be periodic in 3D. This affects the calculation', new_line = .true.) 226 call messages_write('of total energy and forces.') 227 call messages_warning(namespace=namespace) 228 end if 229 230 !%Variable MultiResolutionArea 231 !%Type block 232 !%Section Mesh 233 !%Description 234 !% (Experimental) Multiresolution regions are set with this 235 !% parameter. The first three numbers define the central 236 !% point of the region, and the following ones set 237 !% the radii where resolution changes (measured from the 238 !% central point). 239 !% NOTE: currently, only one area can be set up, and only works in 3D, and in serial. 240 !%End 241 242 if(parse_block(namespace, 'MultiResolutionArea', blk) == 0) then 243 244 call messages_experimental('Multi-resolution') 245 246 if(sb%dim /= 3) call messages_not_implemented('multi-resolution for dim != 3', namespace=namespace) 247 248 ! number of areas 249 sb%hr_area%num_areas = parse_block_n(blk) 250 251 ! number of radii 252 sb%hr_area%num_radii = parse_block_cols(blk, 0) - sb%dim 253 254 sb%hr_area%center = M_ZERO 255 256 ! the central point 257 do idir = 1, sb%dim 258 call parse_block_float(blk, 0, idir - 1, sb%hr_area%center(idir)) 259 end do 260 261 if (sb%hr_area%num_areas /= 1) call messages_input_error(namespace, 'MultiResolutionArea') 262 263 ! the radii 264 SAFE_ALLOCATE(sb%hr_area%radius(1:sb%hr_area%num_radii)) 265 do irad = 1, sb%hr_area%num_radii 266 call parse_block_float(blk, 0, sb%dim + irad - 1, sb%hr_area%radius(irad)) 267 sb%hr_area%radius(irad) = units_to_atomic(units_inp%length, sb%hr_area%radius(irad)) 268 end do 269 call parse_block_end(blk) 270 271 ! Create interpolation points (posi) and weights (ww) 272 273 !%Variable MultiResolutionInterpolationOrder 274 !%Type integer 275 !%Default 5 276 !%Section Mesh 277 !%Description 278 !% The interpolation order in the multiresolution approach (with <tt>MultiResolutionArea</tt>). 279 !%End 280 call messages_obsolete_variable(namespace, 'MR_InterpolationOrder', 'MultiResolutionInterpolationOrder') 281 call parse_variable(namespace, 'MultiResolutionInterpolationOrder', 5, order) 282 call simul_box_interp_init(sb, order, namespace) 283 284 sb%mr_flag = .true. 285 else 286 nullify(sb%hr_area%radius) 287 nullify(sb%hr_area%interp%posi) 288 nullify(sb%hr_area%interp%ww) 289 sb%mr_flag = .false. 290 end if 291 292 POP_SUB(simul_box_init.read_misc) 293 end subroutine read_misc 294 295 296 !-------------------------------------------------------------- 297 subroutine read_box() 298 type(block_t) :: blk 299 300 FLOAT :: default 301 integer :: default_boxshape, idir 302#if defined(HAVE_GDLIB) 303 logical :: found 304 integer :: box_npts 305#endif 306 307 PUSH_SUB(simul_box_init.read_box) 308 ! Read box shape. 309 310 !%Variable BoxShape 311 !%Type integer 312 !%Section Mesh::Simulation Box 313 !%Description 314 !% This variable decides the shape of the simulation box. 315 !% The default is <tt>minimum</tt> for finite systems and <tt>parallelepiped</tt> for periodic systems. 316 !% Note that some incompatibilities apply: 317 !% <ul><li>Spherical or minimum mesh is not allowed for periodic systems. 318 !% <li>Cylindrical mesh is not allowed for systems that are periodic in more than one dimension. 319 !% <li><tt>box_image</tt> is only allowed in 2D.</ul> 320 !%Option sphere 1 321 !% The simulation box will be a sphere of radius <tt>Radius</tt>. (In 2D, this is a circle.) 322 !%Option cylinder 2 323 !% The simulation box will be a cylinder with radius <tt>Radius</tt> and height (in the <i>x</i>-direction) 324 !% of 2 <tt>Xlength</tt>. 325 !%Option minimum 3 326 !% The simulation box will be constructed by adding spheres created around each 327 !% atom (or user-defined potential), of radius <tt>Radius</tt>. 328 !%Option parallelepiped 4 329 !% The simulation box will be a parallelepiped whose dimensions are taken from 330 !% the variable <tt>Lsize</tt>. 331 !%Option box_image 5 332 !% The simulation box will be defined through an image, specified with <tt>BoxShapeImage</tt>. 333 !% White (RGB = 255,255,255) means that the point 334 !% is contained in the simulation box, while any other color means that the point is out. 335 !% The image will be scaled to fit <tt>Lsize</tt>, while its resolution will define the default <tt>Spacing</tt>. 336 !% The actual box may be slightly larger than <tt>Lsize</tt> to ensure one grid point = one pixel for 337 !% default <tt>Spacing</tt>. 338 !%Option user_defined 77 339 !% The shape of the simulation box will be read from the variable <tt>BoxShapeUsDef</tt>. 340 !%Option hypercube 6 341 !% (experimental) The simulation box will be a hypercube or 342 !% hyperparallelepiped. This is equivalent to the 343 !% <tt>parallelepiped</tt> box but it can work with an arbitrary 344 !% number of dimensions. 345 !%End 346 347 if(simul_box_is_periodic(sb)) then 348 default_boxshape = PARALLELEPIPED 349 else 350 default_boxshape = MINIMUM 351 end if 352 call parse_variable(namespace, 'BoxShape', default_boxshape, sb%box_shape) 353 if(.not.varinfo_valid_option('BoxShape', sb%box_shape)) call messages_input_error(namespace, 'BoxShape') 354 select case(sb%box_shape) 355 case(SPHERE, MINIMUM, BOX_USDEF) 356 if(sb%dim > 1 .and. simul_box_is_periodic(sb)) call messages_input_error(namespace, 'BoxShape') 357 case(CYLINDER) 358 if(sb%dim == 2) then 359 message(1) = "BoxShape = cylinder is not meaningful in 2D. Use sphere if you want a circle." 360 call messages_fatal(1, namespace=namespace) 361 end if 362 if(sb%periodic_dim > 1) call messages_input_error(namespace, 'BoxShape') 363 end select 364 365 ! ignore box_shape in 1D 366 if(sb%dim == 1 .and. sb%box_shape /= PARALLELEPIPED .and. sb%box_shape /= HYPERCUBE) & 367 sb%box_shape = SPHERE 368 369 ! Cannot use images in 1D or 3D 370 if(sb%dim /= 2 .and. sb%box_shape == BOX_IMAGE) call messages_input_error(namespace, 'BoxShape') 371 372 if(sb%dim > 3 .and. sb%box_shape /= HYPERCUBE) then 373 message(1) = "For more than 3 dimensions, you can only use the hypercubic box." 374 call messages_fatal(1, namespace=namespace) 375 ! FIXME: why not a hypersphere as another option? 376 ! Also, hypercube should be unified with parallepiped. 377 end if 378 379 sb%rsize = -M_ONE 380 !%Variable Radius 381 !%Type float 382 !%Section Mesh::Simulation Box 383 !%Description 384 !% Defines the radius for <tt>BoxShape</tt> = <tt>sphere</tt>, 385 !% <tt>cylinder</tt>, or <tt>minimum</tt>. Must be a positive 386 !% number. If not specified, the code will look for values in 387 !% the <tt>Species</tt> block, or, from the default 388 !% pseudopotential parameters. In these cases, for 389 !% <tt>minimum</tt>, a different radius is used for each 390 !% species, while for other shapes, the maximum radius is used. 391 !%End 392 select case(sb%box_shape) 393 case(SPHERE, CYLINDER) 394 call parse_variable(namespace, 'Radius', def_rsize, sb%rsize, units_inp%length) 395 if(sb%rsize < M_ZERO) call messages_input_error(namespace, 'radius') 396 if(def_rsize>M_ZERO) call messages_check_def(sb%rsize, .false., def_rsize, 'radius', units_out%length) 397 case(MINIMUM) 398 399 if(geo%reduced_coordinates) then 400 message(1) = "The 'minimum' box shape cannot be used if atomic positions" 401 message(2) = "are given as reduced coordinates." 402 call messages_fatal(2, namespace=namespace) 403 end if 404 405 default=sb%rsize 406 call parse_variable(namespace, 'radius', default, sb%rsize, units_inp%length) 407 if(sb%rsize < M_ZERO .and. def_rsize < M_ZERO) call messages_input_error(namespace, 'Radius') 408 end select 409 410 if(sb%box_shape == CYLINDER) then 411 !%Variable Xlength 412 !%Default <tt>Radius</tt> 413 !%Type float 414 !%Section Mesh::Simulation Box 415 !%Description 416 !% If <tt>BoxShape</tt> is <tt>cylinder</tt>, the total length of the cylinder is twice <tt>Xlength</tt>. 417 !%End 418 if(sb%rsize > M_ZERO) then 419 default = sb%rsize 420 else 421 default = def_rsize 422 end if 423 424 call parse_variable(namespace, 'Xlength', default, sb%xsize, units_inp%length) 425 if(def_rsize > M_ZERO .and. sb%periodic_dim == 0) & 426 call messages_check_def(sb%xsize, .false., def_rsize, 'xlength', units_out%length) 427 end if 428 429 sb%lsize = M_ZERO 430 if(sb%box_shape == PARALLELEPIPED .or. sb%box_shape == HYPERCUBE .or. & 431 sb%box_shape == BOX_IMAGE .or. sb%box_shape == BOX_USDEF) then 432 433 !%Variable Lsize 434 !%Type block 435 !%Section Mesh::Simulation Box 436 !%Description 437 !% If <tt>BoxShape</tt> is <tt>parallelepiped</tt>, <tt>hypercube</tt>, 438 !% <tt>box_image</tt>, or <tt>user_defined</tt>, this is a 439 !% block of the form: 440 !% 441 !% <tt>%Lsize 442 !% <br> sizex | sizey | sizez | ... 443 !% <br>%</tt> 444 !% 445 !% where the <tt>size*</tt> are half the lengths of the box in each direction. 446 !% 447 !% The number of columns must match the dimensionality of the 448 !% calculation. If you want a cube you can also set <tt>Lsize</tt> as a 449 !% single variable. 450 !%End 451 452 if(all(geo%lsize(1:sb%dim) > M_ZERO)) then 453 ! use value read from XSF lattice vectors 454 sb%lsize(:) = geo%lsize(:) 455 else if(parse_block(namespace, 'Lsize', blk) == 0) then 456 if(parse_block_cols(blk,0) < sb%dim .and. .not. parse_is_defined(namespace, 'LatticeVectors')) & 457 call messages_input_error(namespace, 'Lsize') 458 do idir = 1, sb%dim 459 call parse_block_float(blk, 0, idir - 1, sb%lsize(idir), units_inp%length) 460 if(def_rsize > M_ZERO .and. sb%periodic_dim < idir) & 461 call messages_check_def(sb%lsize(idir), .false., def_rsize, 'Lsize', units_out%length) 462 end do 463 call parse_block_end(blk) 464 else if ((parse_is_defined(namespace, 'Lsize'))) then 465 call parse_variable(namespace, 'Lsize', -M_ONE, sb%lsize(1), units_inp%length) 466 if(abs(sb%lsize(1)+M_ONE) <= M_EPSILON) then 467 call messages_input_error(namespace, 'Lsize') 468 end if 469 if(def_rsize > M_ZERO .and. sb%periodic_dim < sb%dim) & 470 call messages_check_def(sb%lsize(1), .false., def_rsize, 'Lsize', units_out%length) 471 sb%lsize(1:sb%dim) = sb%lsize(1) 472 else 473 message(1) = "Lsize was not found in input file. Continuing anyway." 474 call messages_warning(1, namespace=namespace) 475 end if 476 else 477 ! if not a compatible box-shape 478 if(all(geo%lsize(1:sb%dim) > M_ZERO)) then 479 message(1) = "Ignoring lattice vectors from XSF file." 480 call messages_warning(1, namespace=namespace) 481 end if 482 end if 483 484 ! read in image for box_image 485 if(sb%box_shape == BOX_IMAGE) then 486 487 !%Variable BoxShapeImage 488 !%Type string 489 !%Section Mesh::Simulation Box 490 !%Description 491 !% Name of the file that contains the image that defines the simulation box 492 !% when <tt>BoxShape = box_image</tt>. No default. Will search in current 493 !% directory and <tt>OCTOPUS-HOME/share/</tt>. 494 !%End 495#if defined(HAVE_GDLIB) 496 call parse_variable(namespace, 'BoxShapeImage', '', sb%filename) 497 if(trim(sb%filename) == "") then 498 message(1) = "Must specify BoxShapeImage if BoxShape = box_image." 499 call messages_fatal(1, namespace=namespace) 500 end if 501 502 ! Find out the file and read it. 503 inquire(file=trim(sb%filename), exist=found) 504 if(.not. found) then 505 message(1) = "Could not find file '" // trim(sb%filename) // "' for BoxShape = box_image." 506 507 sb%filename = trim(conf%share) // '/' // trim(sb%filename) 508 inquire(file=trim(sb%filename), exist=found) 509 510 if(.not. found) call messages_fatal(1, namespace=namespace) 511 end if 512 513 sb%image = gdlib_image_create_from(sb%filename) 514 if(.not.c_associated(sb%image)) then 515 message(1) = "Could not open file '" // trim(sb%filename) // "' for BoxShape = box_image." 516 call messages_fatal(1, namespace=namespace) 517 end if 518 sb%image_size(1) = gdlib_image_sx(sb%image) 519 sb%image_size(2) = gdlib_image_sy(sb%image) 520 521 ! adjust Lsize if necessary to ensure that one grid point = one pixel 522 do idir = 1, 2 523 box_npts = sb%image_size(idir) 524 if((idir > sb%periodic_dim .and. even(sb%image_size(idir))) .or. & 525 (idir <= sb%periodic_dim .and. odd(sb%image_size(idir)))) then 526 box_npts = box_npts + 1 527 sb%lsize(idir) = sb%lsize(idir) * box_npts / sb%image_size(idir) 528 end if 529 end do 530#else 531 message(1) = "To use 'BoxShape = box_image', you have to compile Octopus" 532 message(2) = "with GD library support." 533 call messages_fatal(2, namespace=namespace) 534#endif 535 end if 536 537 ! read in box shape for user-defined boxes 538 if(sb%box_shape == BOX_USDEF) then 539 540 !%Variable BoxShapeUsDef 541 !%Type string 542 !%Section Mesh::Simulation Box 543 !%Description 544 !% Boolean expression that defines the interior of the simulation box. For example, 545 !% <tt>BoxShapeUsDef = "(sqrt(x^2+y^2) <= 4) && z>-2 && z<2"</tt> defines a cylinder 546 !% with axis parallel to the <i>z</i>-axis. 547 !%End 548 549 call parse_variable(namespace, 'BoxShapeUsDef', 'x^2+y^2+z^2 < 4', sb%user_def) 550 call conv_to_C_string(sb%user_def) 551 end if 552 553 ! fill in lsize structure 554 select case(sb%box_shape) 555 case(SPHERE) 556 sb%lsize(1:sb%dim) = sb%rsize 557 case(CYLINDER) 558 sb%lsize(1) = sb%xsize 559 sb%lsize(2:sb%dim) = sb%rsize 560 case(MINIMUM) 561 do idir = 1, sb%dim 562 if(sb%rsize > M_ZERO) then 563 sb%lsize(idir) = maxval(abs(geo%atom(:)%x(idir))) + sb%rsize 564 else 565 sb%lsize(idir) = maxval(abs(geo%atom(:)%x(idir))) + def_rsize 566 end if 567 end do 568 end select 569 570 call messages_obsolete_variable(namespace, 'BoxOffset') 571 572 POP_SUB(simul_box_init.read_box) 573 end subroutine read_box 574 575 end subroutine simul_box_init 576 577 ! ------------------------------------------------------------ 578 subroutine simul_box_lookup_init(this, geo) 579 type(simul_box_t), intent(inout) :: this 580 type(geometry_t), intent(in) :: geo 581 ! 582 FLOAT, allocatable :: pos(:, :) 583 integer :: iatom 584 ! 585 PUSH_SUB(simul_box_lookup_init) 586 587 SAFE_ALLOCATE(pos(1:this%dim,1:geo%natoms)) 588 589 do iatom = 1, geo%natoms 590 pos(:,iatom) = geo%atom(iatom)%x(1:this%dim) 591 end do 592 593 call lookup_init(this%atom_lookup, this%dim, geo%natoms, pos) 594 595 SAFE_DEALLOCATE_A(pos) 596 POP_SUB(simul_box_lookup_init) 597 return 598 end subroutine simul_box_lookup_init 599 600 ! ------------------------------------------------------------ 601 subroutine simul_box_interp_init(this, order, namespace) 602 type(simul_box_t), intent(inout) :: this 603 integer, intent(in) :: order 604 type(namespace_t), intent(in) :: namespace 605 ! 606 FLOAT, allocatable, dimension(:) :: pos 607 integer :: ii 608 ! 609 PUSH_SUB(simul_box_interp_init) 610 this%hr_area%interp%order=order 611 if(this%hr_area%interp%order<=0) then 612 message(1) = "The value for MultiResolutionInterpolationOrder must be > 0." 613 call messages_fatal(1, namespace=namespace) 614 end if 615 this%hr_area%interp%nn=2*this%hr_area%interp%order 616 SAFE_ALLOCATE(pos(1:this%hr_area%interp%nn)) 617 SAFE_ALLOCATE(this%hr_area%interp%ww(1:this%hr_area%interp%nn)) 618 SAFE_ALLOCATE(this%hr_area%interp%posi(1:this%hr_area%interp%nn)) 619 do ii = 1, this%hr_area%interp%order 620 this%hr_area%interp%posi(ii)=1+2*(ii-1) 621 this%hr_area%interp%posi(this%hr_area%interp%order+ii)=-this%hr_area%interp%posi(ii) 622 pos(ii)=this%hr_area%interp%posi(ii) 623 pos(this%hr_area%interp%order+ii)=-pos(ii) 624 end do 625 call interpolation_coefficients(this%hr_area%interp%nn, pos, M_ZERO, this%hr_area%interp%ww) 626 SAFE_DEALLOCATE_A(pos) 627 POP_SUB(simul_box_interp_init) 628 return 629 end subroutine simul_box_interp_init 630 631 !-------------------------------------------------------------- 632 subroutine simul_box_build_lattice(sb, namespace, rlattice_primitive) 633 type(simul_box_t), intent(inout) :: sb 634 type(namespace_t), intent(in) :: namespace 635 FLOAT, optional, intent(in) :: rlattice_primitive(:,:) 636 637 type(block_t) :: blk 638 FLOAT :: norm, lparams(3) 639 integer :: idim, jdim, ncols 640 logical :: has_angles 641 FLOAT :: angles(1:MAX_DIM), cosang, a2, aa, cc 642 FLOAT, parameter :: tol_angle = CNST(1.0e-6) 643 644 PUSH_SUB(simul_box_build_lattice) 645 646 sb%alpha = CNST(90.0) 647 sb%beta = CNST(90.0) 648 sb%gamma = CNST(90.0) 649 650 if(present(rlattice_primitive)) then 651 sb%rlattice_primitive(1:sb%dim, 1:sb%dim) = rlattice_primitive(1:sb%dim, 1:sb%dim) 652 sb%nonorthogonal = .false. 653 else 654 655 656 !%Variable LatticeParameters 657 !%Type block 658 !%Default 1 | 1 | 1 659 !%Section Mesh::Simulation Box 660 !%Description 661 !% The lattice parameters (a, b, c). 662 !% This option is incompatible with Lsize and either one of the 663 !% two must be specified in the input file for periodic systems. 664 !% A second optional line can be used tu define the angles between the lattice vectors 665 !%End 666 lparams(:) = M_ONE 667 has_angles = .false. 668 angles = CNST(90.0) 669 670 if (parse_block(namespace, 'LatticeParameters', blk) == 0) then 671 do idim = 1, sb%dim 672 call parse_block_float(blk, 0, idim - 1, lparams(idim)) 673 end do 674 675 if(parse_block_n(blk) > 1) then ! we have a shift, or even more 676 ncols = parse_block_cols(blk, 1) 677 if(ncols /= sb%dim) then 678 write(message(1),'(a,i3,a,i3)') 'LatticeParameters angle has ', ncols, ' columns but must have ', sb%dim 679 call messages_fatal(1, namespace=namespace) 680 end if 681 do idim = 1, sb%dim 682 call parse_block_float(blk, 1, idim - 1, angles(idim)) 683 end do 684 has_angles = .true. 685 end if 686 call parse_block_end(blk) 687 688 if (parse_is_defined(namespace, 'Lsize')) then 689 message(1) = 'LatticeParameters is incompatible with Lsize' 690 call messages_print_var_info(stdout, "LatticeParameters") 691 call messages_fatal(1, namespace=namespace) 692 end if 693 694 end if 695 696 if( has_angles ) then 697 sb%alpha = angles(1) 698 sb%beta = angles(2) 699 sb%gamma = angles(3) 700 701 !Converting the angles to LatticeVectors 702 !See 57_iovars/ingeo.F90 in Abinit for details 703 if( abs(angles(1)-angles(2))< tol_angle .and. abs(angles(2)-angles(3))< tol_angle .and. & 704 (abs(angles(1)-CNST(90.0))+abs(angles(2)-CNST(90.0))+abs(angles(3)-CNST(90.0)))> tol_angle ) then 705 706 cosang=cos(M_PI*angles(1)/CNST(180.0)); 707 a2=M_TWO/M_THREE*(M_ONE-cosang); 708 aa=sqrt(a2); 709 cc=sqrt(M_ONE-a2); 710 sb%rlattice_primitive(1,1) = aa 711 sb%rlattice_primitive(2,1) = M_ZERO 712 sb%rlattice_primitive(3,1) = cc 713 sb%rlattice_primitive(1,2) =-M_HALF*aa 714 sb%rlattice_primitive(2,2) = M_HALF*sqrt(M_THREE)*aa 715 sb%rlattice_primitive(3,2) = cc 716 sb%rlattice_primitive(1,3) =-M_HALF*aa 717 sb%rlattice_primitive(2,3) =-M_HALF*sqrt(M_THREE)*aa 718 sb%rlattice_primitive(3,3) = cc 719 else 720 sb%rlattice_primitive(1,1) = M_ONE 721 sb%rlattice_primitive(2,1) = M_ZERO 722 sb%rlattice_primitive(3,1) = M_ZERO 723 sb%rlattice_primitive(1,2) = cos(M_PI*angles(3)/CNST(180.0)) 724 sb%rlattice_primitive(2,2) = sin(M_PI*angles(3)/CNST(180.0)) 725 sb%rlattice_primitive(3,2) = M_ZERO 726 sb%rlattice_primitive(1,3) = cos(M_PI*angles(2)/CNST(180.0)) 727 sb%rlattice_primitive(2,3) = (cos(M_PI*angles(1)/CNST(180.0))-sb%rlattice_primitive(1,2)* sb%rlattice_primitive(1,3))& 728 /sb%rlattice_primitive(2,2) 729 sb%rlattice_primitive(3,3) = sqrt(M_ONE-sb%rlattice_primitive(1,3)**2-sb%rlattice_primitive(2,3)**2) 730 end if 731 732 if (parse_is_defined(namespace, 'LatticeVectors')) then 733 message(1) = 'LatticeParameters with angles is incompatible with LatticeVectors' 734 call messages_print_var_info(stdout, "LatticeParameters") 735 call messages_fatal(1, namespace=namespace) 736 end if 737 738 if(any(abs(angles-CNST(90.0)) > M_EPSILON )) then 739 sb%nonorthogonal = .true. 740 else 741 sb%nonorthogonal = .false. 742 end if 743 744 else 745 746 !%Variable LatticeVectors 747 !%Type block 748 !%Default simple cubic 749 !%Section Mesh::Simulation Box 750 !%Description 751 !% Primitive lattice vectors. Vectors are stored in rows. 752 !% Default: 753 !% <br><br><tt>%LatticeVectors 754 !% <br> 1.0 | 0.0 | 0.0 755 !% <br> 0.0 | 1.0 | 0.0 756 !% <br> 0.0 | 0.0 | 1.0 757 !% <br>%<br></tt> 758 !%End 759 sb%rlattice_primitive = M_ZERO 760 sb%nonorthogonal = .false. 761 do idim = 1, sb%dim 762 sb%rlattice_primitive(idim, idim) = M_ONE 763 end do 764 765 if (parse_block(namespace, 'LatticeVectors', blk) == 0) then 766 do idim = 1, sb%dim 767 do jdim = 1, sb%dim 768 call parse_block_float(blk, idim - 1, jdim - 1, sb%rlattice_primitive(jdim, idim)) 769 if(idim /= jdim .and. abs(sb%rlattice_primitive(jdim, idim)) > M_EPSILON) sb%nonorthogonal = .true. 770 enddo 771 end do 772 call parse_block_end(blk) 773 774 end if 775 end if 776 777 ! Always need Lsize for periodic systems even if LatticeVectors block is not present 778 if (.not. parse_is_defined(namespace, 'Lsize') .and. sb%periodic_dim > 0) then 779 do idim = 1, sb%dim 780 if (sb%lsize(idim) == M_ZERO) then 781 sb%lsize(idim) = lparams(idim)*M_HALF 782 end if 783 end do 784 end if 785 786 end if 787 788 sb%rlattice = M_ZERO 789 do idim = 1, sb%dim 790 norm = sqrt(sum(sb%rlattice_primitive(1:sb%dim, idim)**2)) 791 sb%lsize(idim) = sb%lsize(idim) * norm 792 do jdim = 1, sb%dim 793 sb%rlattice_primitive(jdim, idim) = sb%rlattice_primitive(jdim, idim) / norm 794 sb%rlattice(jdim, idim) = sb%rlattice_primitive(jdim, idim) * M_TWO*sb%lsize(idim) 795 end do 796 end do 797 798 call reciprocal_lattice(sb%rlattice, sb%klattice, sb%rcell_volume, sb%dim, namespace) 799 sb%klattice = sb%klattice * M_TWO*M_PI 800 801 call reciprocal_lattice(sb%rlattice_primitive, sb%klattice_primitive, sb%volume_element, sb%dim, namespace) 802 803 if(sb%dim == 3) then 804 sb%surface_element(1) = sqrt(abs(sum(dcross_product(sb%rlattice_primitive(1:3, 2), sb%rlattice_primitive(1:3, 3))**2))) 805 sb%surface_element(2) = sqrt(abs(sum(dcross_product(sb%rlattice_primitive(1:3, 3), sb%rlattice_primitive(1:3, 1))**2))) 806 sb%surface_element(3) = sqrt(abs(sum(dcross_product(sb%rlattice_primitive(1:3, 1), sb%rlattice_primitive(1:3, 2))**2))) 807 end if 808 809 ! rlattice_primitive is the A matrix from Chelikowski PRB 78 075109 (2008) 810 ! klattice_primitive is the transpose (!) of the B matrix, with no 2 pi factor included 811 ! klattice is the proper reciprocal lattice vectors, with 2 pi factor, and in units of 1/bohr 812 ! The F matrix of Chelikowski is matmul(transpose(sb%klattice_primitive), sb%klattice_primitive) 813 sb%rmetric = matmul(transpose(sb%rlattice_primitive), sb%rlattice_primitive) 814 if(.not. has_angles .and. sb%dim == 3) then 815 !We compute the angles from the lattice vectors 816 sb%alpha=acos(sb%rmetric(2,3)/sqrt(sb%rmetric(2,2)*sb%rmetric(3,3)))/M_PI*CNST(180.0) 817 sb%beta =acos(sb%rmetric(1,3)/sqrt(sb%rmetric(1,1)*sb%rmetric(3,3)))/M_PI*CNST(180.0) 818 sb%gamma=acos(sb%rmetric(1,2)/sqrt(sb%rmetric(1,1)*sb%rmetric(2,2)))/M_PI*CNST(180.0) 819 end if 820 821 POP_SUB(simul_box_build_lattice) 822 end subroutine simul_box_build_lattice 823 824 825 !> This function adjusts the coordinates defined in the geometry 826 !! object. If coordinates were given in reduced coordinates it 827 !! converts them to real coordinates and it checks that the atoms 828 !! are inside the box. 829 !! 830 !! If atoms are not in the box: if the system is periodic, the atoms 831 !! are moved inside the box, if the system is finite, nothing 832 !! happens or a warning is written, depending on the argument 833 !! warn_if_not. 834 ! --------------------------------------------------------- 835 subroutine simul_box_atoms_in_box(sb, geo, namespace, warn_if_not, die_if_not) 836 type(simul_box_t), intent(in) :: sb 837 type(geometry_t), intent(inout) :: geo 838 type(namespace_t), intent(in) :: namespace 839 logical, intent(in) :: warn_if_not 840 logical, optional, intent(in) :: die_if_not 841 842 integer :: iatom, pd 843 logical :: die_if_not_ 844 845 PUSH_SUB(simul_box_atoms_in_box) 846 847 die_if_not_ = optional_default(die_if_not, .false.) 848 pd = sb%periodic_dim 849 850 do iatom = 1, geo%natoms 851 852 call simul_box_periodic_atom_in_box(sb, geo, geo%atom(iatom)%x(:)) 853 854 if(geo%reduced_coordinates) then 855 geo%atom(iatom)%x(pd + 1:sb%dim) = M_TWO*sb%lsize(pd + 1:sb%dim)*geo%atom(iatom)%x(pd + 1:sb%dim) 856 end if 857 858 if( .not. simul_box_in_box(sb, geo, geo%atom(iatom)%x, namespace) ) then 859 write(message(1), '(a,i5,a)') "Atom ", iatom, " is outside the box." 860 if (sb%periodic_dim /= sb%dim) then 861 ! FIXME: This could fail for partial periodicity systems 862 ! because simul_box_in_box is too strict with atoms close to 863 ! the upper boundary to the cell. 864 if(warn_if_not) call messages_warning(1, namespace=namespace) 865 if(die_if_not_) call messages_fatal(1, namespace=namespace) 866 end if 867 end if 868 869 end do 870 871 ! done with the conversion to real coordinates 872 geo%reduced_coordinates = .false. 873 874 POP_SUB(simul_box_atoms_in_box) 875 end subroutine simul_box_atoms_in_box 876 877 ! -------------------------------------------------------- 878 879 subroutine simul_box_periodic_atom_in_box(sb, geo, ratom) 880 type(simul_box_t), intent(in) :: sb 881 type(geometry_t), intent(in) :: geo 882 FLOAT, intent(inout) :: ratom(:) 883 884 FLOAT :: xx(1:MAX_DIM) 885 integer :: pd, idir 886 887 pd = sb%periodic_dim 888 889 if (simul_box_is_periodic(sb)) then 890 if(.not. geo%reduced_coordinates) then 891 !convert the position to reduced coordinates 892 xx(1:pd) = matmul(ratom(1:pd), sb%klattice(1:pd, 1:pd))/(M_TWO*M_PI) 893 else 894 ! in this case coordinates are already in reduced space 895 xx(1:pd) = ratom(1:pd) 896 end if 897 898 xx(1:pd) = xx(1:pd) + M_HALF 899 do idir = 1, pd 900 xx(idir) = xx(idir) - anint(xx(idir)) 901 if(xx(idir) < -CNST(1.0e-6)) & 902 xx(idir) = xx(idir) + M_ONE 903 end do 904 ASSERT(all(xx(1:pd) >= -CNST(1.0e-6))) 905 ASSERT(all(xx(1:pd) < CNST(1.0))) 906 907 xx(1:pd) = (xx(1:pd) - M_HALF) 908 ratom(1:pd) = matmul(sb%rlattice(1:pd, 1:pd), xx(1:pd)) 909 910 end if 911 912 913 end subroutine simul_box_periodic_atom_in_box 914 915 !-------------------------------------------------------------- 916 subroutine reciprocal_lattice(rv, kv, volume, dim, namespace) 917 FLOAT, intent(in) :: rv(:,:) !< (1:MAX_DIM, 1:MAX_DIM) 918 FLOAT, intent(out) :: kv(:,:) !< (1:MAX_DIM, 1:MAX_DIM) 919 FLOAT, intent(out) :: volume 920 integer, intent(in) :: dim 921 type(namespace_t), intent(in) :: namespace 922 923 integer :: ii 924 FLOAT :: cross(1:3), rv3(1:3, 1:3) 925 926 PUSH_SUB(reciprocal_lattice) 927 928 kv(:,:) = M_ZERO 929 930 select case(dim) 931 case(3) 932 cross(1:3) = dcross_product(rv(1:3, 2), rv(1:3, 3)) 933 volume = dot_product(rv(1:3, 1), cross(1:3)) 934 935 kv(1:3, 1) = dcross_product(rv(:, 2), rv(:, 3))/volume 936 kv(1:3, 2) = dcross_product(rv(:, 3), rv(:, 1))/volume 937 kv(1:3, 3) = dcross_product(rv(:, 1), rv(:, 2))/volume 938 case(2) 939 rv3(1:3, 1:3) = M_ZERO 940 rv3(1:2, 1:2) = rv(1:2, 1:2) 941 rv3(3, 3) = M_ONE 942 cross(1:3) = dcross_product(rv3(1:3, 1), rv3(1:3, 2)) 943 volume = dot_product(rv3(1:3, 3), cross(1:3)) 944 945 kv(1:3, 1) = dcross_product(rv3(:, 2), rv3(:, 3))/volume 946 kv(1:3, 2) = dcross_product(rv3(:, 3), rv3(:, 1))/volume 947 case(1) 948 volume = rv(1, 1) 949 kv(1, 1) = M_ONE / rv(1, 1) 950 case default ! dim > 3 951 message(1) = "Reciprocal lattice for dim > 3 assumes no periodicity." 952 call messages_warning(1, namespace=namespace) 953 volume = M_ONE 954 do ii = 1, dim 955 kv(ii, ii) = M_ONE/rv(ii,ii) 956 ! At least initialize the thing 957 volume = volume * sqrt(sum(rv(:, ii)**2)) 958 end do 959 end select 960 961 if ( volume < M_ZERO ) then 962 message(1) = "Your lattice vectors form a left-handed system." 963 call messages_fatal(1, namespace=namespace) 964 end if 965 966 POP_SUB(reciprocal_lattice) 967 end subroutine reciprocal_lattice 968 969 !-------------------------------------------------------------- 970 subroutine simul_box_end(sb) 971 type(simul_box_t), intent(inout) :: sb 972 973 PUSH_SUB(simul_box_end) 974 975 call symmetries_end(sb%symm) 976 977 call lookup_end(sb%atom_lookup) 978 call kpoints_end(sb%kpoints) 979 980 SAFE_DEALLOCATE_P(sb%hr_area%radius) 981 SAFE_DEALLOCATE_P(sb%hr_area%interp%ww) 982 SAFE_DEALLOCATE_P(sb%hr_area%interp%posi) 983 984#ifdef HAVE_GDLIB 985 if(sb%box_shape == BOX_IMAGE) & 986 call gdlib_imagedestroy(sb%image) 987#endif 988 989 POP_SUB(simul_box_end) 990 end subroutine simul_box_end 991 992 993 !-------------------------------------------------------------- 994 recursive subroutine simul_box_write_info(sb, geo, iunit) 995 type(simul_box_t), intent(in) :: sb 996 type(geometry_t), intent(in) :: geo 997 integer, intent(in) :: iunit 998 999 character(len=15), parameter :: bs(6) = (/ & 1000 'sphere ', & 1001 'cylinder ', & 1002 'minimum ', & 1003 'parallelepiped', & 1004 'image-defined ', & 1005 'hypercube '/) 1006 1007 integer :: idir, idir2, ispec 1008 1009 PUSH_SUB(simul_box_write_info) 1010 1011 write(message(1),'(a)') 'Simulation Box:' 1012 if(sb%box_shape == BOX_USDEF) then 1013 write(message(2), '(a)') ' Type = user-defined' 1014 else if(sb%box_shape == BOX_IMAGE) then 1015 write(message(2), '(3a,i6,a,i6)') ' Type = defined by image "', trim(sb%filename), '", ', & 1016 sb%image_size(1), ' x ', sb%image_size(2) 1017 else 1018 write(message(2), '(2a)') ' Type = ', trim(bs(sb%box_shape)) 1019 end if 1020 call messages_info(2, iunit) 1021 1022 if(sb%box_shape == SPHERE .or. sb%box_shape == CYLINDER & 1023 .or. (sb%box_shape == MINIMUM .and. sb%rsize > M_ZERO)) then 1024 write(message(1), '(3a,f7.3)') ' Radius [', trim(units_abbrev(units_out%length)), '] = ', & 1025 units_from_atomic(units_out%length, sb%rsize) 1026 call messages_info(1, iunit) 1027 end if 1028 1029 if (sb%box_shape == MINIMUM .and. sb%rsize <= M_ZERO) then 1030 do ispec = 1, geo%nspecies 1031 write(message(1), '(a,a5,5x,a,f7.3,2a)') ' Species = ', trim(species_label(geo%species(ispec))), 'Radius = ', & 1032 units_from_atomic(units_out%length, species_def_rsize(geo%species(ispec))), ' ', trim(units_abbrev(units_out%length)) 1033 call messages_info(1, iunit) 1034 end do 1035 end if 1036 1037 if(sb%box_shape == CYLINDER) then 1038 write(message(1), '(3a,f7.3)') ' Xlength [', trim(units_abbrev(units_out%length)), '] = ', & 1039 units_from_atomic(units_out%length, sb%xsize) 1040 call messages_info(1, iunit) 1041 end if 1042 1043 if(sb%box_shape == PARALLELEPIPED) then 1044 write(message(1),'(3a, 99(a, f8.3), a)') & 1045 ' Lengths [', trim(units_abbrev(units_out%length)), '] = ', & 1046 '(', (units_from_atomic(units_out%length, M_TWO*sb%lsize(idir)), ',', idir = 1, sb%dim - 1), & 1047 units_from_atomic(units_out%length, M_TWO*sb%lsize(sb%dim)), ')' 1048 call messages_info(1, iunit) 1049 end if 1050 1051 write(message(1), '(a,i1,a)') ' Octopus will run in ', sb%dim, ' dimension(s).' 1052 write(message(2), '(a,i1,a)') ' Octopus will treat the system as periodic in ', & 1053 sb%periodic_dim, ' dimension(s).' 1054 call messages_info(2, iunit) 1055 1056 if(sb%periodic_dim > 0 .or. sb%box_shape == PARALLELEPIPED) then 1057 write(message(1),'(1x)') 1058 write(message(2),'(a,3a,a)') ' Lattice Vectors [', trim(units_abbrev(units_out%length)), ']' 1059 do idir = 1, sb%dim 1060 write(message(2+idir),'(9f12.6)') (units_from_atomic(units_out%length, sb%rlattice(idir2, idir)), & 1061 idir2 = 1, sb%dim) 1062 end do 1063 call messages_info(2+sb%dim, iunit) 1064 1065 write(message(1),'(a,f18.4,3a,i1.1,a)') & 1066 ' Cell volume = ', units_from_atomic(units_out%length**sb%dim, sb%rcell_volume), & 1067 ' [', trim(units_abbrev(units_out%length**sb%dim)), ']' 1068 call messages_info(1, iunit) 1069 1070 write(message(1),'(a,3a,a)') ' Reciprocal-Lattice Vectors [', trim(units_abbrev(units_out%length**(-1))), ']' 1071 do idir = 1, sb%dim 1072 write(message(1+idir),'(3f12.6)') (units_from_atomic(unit_one / units_out%length, sb%klattice(idir2, idir)), & 1073 idir2 = 1, sb%dim) 1074 end do 1075 call messages_info(1+sb%dim, iunit) 1076 1077 if(sb%dim == 3) then 1078 write(message(1),'(a)') ' Cell angles [degree]' 1079 write(message(2),'(a, f8.3)') ' alpha = ', sb%alpha 1080 write(message(3),'(a, f8.3)') ' beta = ', sb%beta 1081 write(message(4),'(a, f8.3)') ' gamma = ', sb%gamma 1082 call messages_info(4, iunit) 1083 end if 1084 end if 1085 1086 POP_SUB(simul_box_write_info) 1087 end subroutine simul_box_write_info 1088 1089 subroutine simul_box_write_short_info(sb, iunit) 1090 type(simul_box_t), intent(in) :: sb 1091 integer, intent(in) :: iunit 1092 1093 integer :: idir1, idir2 1094 1095 PUSH_SUB(simul_box_write_short_info) 1096 1097 write(iunit, '(a,i1,a)', advance='no') 'Dimensions = ', sb%dim, '; ' 1098 write(iunit, '(a,i1,a)', advance='no') 'PeriodicDimensions = ', sb%periodic_dim, '; ' 1099 1100 write(iunit, '(a)', advance='no') 'BoxShape = ' 1101 select case(sb%box_shape) 1102 case(SPHERE) 1103 write(iunit, '(a,f11.6,a)', advance='no') 'sphere; Radius =', units_from_atomic(unit_angstrom, sb%rsize), ' Ang' 1104 case(CYLINDER) 1105 write(iunit, '(a,f11.6,a,f11.6,a)', advance='no') 'cylinder, Radius =', units_from_atomic(unit_angstrom, sb%rsize), & 1106 ' Ang; Xlength =', units_from_atomic(unit_angstrom, sb%xsize), ' Ang' 1107 case(MINIMUM) 1108 write(iunit, '(a,f11.6,a)', advance='no') 'minimum; Radius =', units_from_atomic(unit_angstrom, sb%rsize), ' Ang' 1109 case(PARALLELEPIPED) 1110 write(iunit, '(a)', advance='no') 'parallelepiped; LatticeVectors[Ang] = ' 1111 do idir1 = 1, sb%dim 1112 write(iunit, '(a)', advance='no') '[' 1113 do idir2 = 1, sb%dim 1114 write(iunit, '(x,f11.6)', advance='no') units_from_atomic(unit_angstrom, sb%rlattice(idir2, idir1)) 1115 end do 1116 write(iunit, '(a)', advance='no') ']' 1117 if(idir1 < sb%dim) then 1118 write(iunit, '(a)', advance='no') ', ' 1119 end if 1120 end do 1121 case(BOX_IMAGE) 1122 write(iunit, '(a)', advance='no') 'box_image; BoxShapeImage = '//trim(sb%filename) 1123 case(HYPERCUBE) 1124 write(iunit, '(a)', advance='no') 'hypercube' ! add parameters? 1125 case(BOX_USDEF) 1126 write(iunit, '(a)', advance='no') 'user_defined; BoxShapeUsDef = "'//trim(sb%user_def)//'"' 1127 end select 1128 1129 write(iunit, '()') 1130 POP_SUB(simul_box_write_short_info) 1131 1132 end subroutine simul_box_write_short_info 1133 1134 !-------------------------------------------------------------- 1135 !> Checks if a mesh point belongs to the actual mesh. 1136 logical function simul_box_in_box(sb, geo, yy, namespace) result(in_box) 1137 type(simul_box_t), intent(in) :: sb 1138 type(geometry_t), intent(in) :: geo 1139 FLOAT, intent(in) :: yy(:) 1140 type(namespace_t), intent(in) :: namespace 1141 1142 FLOAT :: xx(1:MAX_DIM, 1) 1143 logical :: in_box2(1) 1144 1145 ! no push_sub because this function is called very frequently 1146 1147 xx(1:sb%dim, 1) = yy(1:sb%dim) 1148 1149 call simul_box_in_box_vec(sb, geo, 1, xx, in_box2, namespace) 1150 in_box = in_box2(1) 1151 1152 end function simul_box_in_box 1153 1154 1155 !-------------------------------------------------------------- 1156 !> Checks if a group of mesh points belong to the actual mesh. 1157 subroutine simul_box_in_box_vec(sb, geo, npoints, point, in_box, namespace) 1158 type(simul_box_t), intent(in) :: sb 1159 type(geometry_t), intent(in) :: geo 1160 integer, intent(in) :: npoints 1161 FLOAT, intent(in) :: point(:, :) 1162 logical, intent(out) :: in_box(:) 1163 type(namespace_t), intent(in) :: namespace 1164 1165 FLOAT, parameter :: DELTA = CNST(1e-12) 1166 FLOAT :: rr, re, im, dist2, radius 1167 FLOAT :: llimit(MAX_DIM), ulimit(MAX_DIM) 1168 FLOAT, allocatable :: xx(:, :) 1169 integer :: ip, idir, iatom, ilist 1170 integer, allocatable :: nlist(:) 1171 integer, pointer :: list(:, :) 1172 1173#if defined(HAVE_GDLIB) 1174 integer :: red, green, blue, ix, iy 1175#endif 1176 1177 ! no push_sub because this function is called very frequently 1178 SAFE_ALLOCATE(xx(1:sb%dim, 1:npoints)) 1179 xx = M_ZERO 1180 1181 !convert from Cartesian to reduced lattice coord 1182 if(npoints == 1) then 1183 xx(1:sb%dim, 1) = matmul(point(1:sb%dim, 1), sb%klattice_primitive(1:sb%dim, 1:sb%dim)) 1184 else 1185 call lalg_gemmt(sb%dim, npoints, sb%dim, M_ONE, sb%klattice_primitive, point, M_ZERO, xx) 1186 end if 1187 1188 select case(sb%box_shape) 1189 case(SPHERE) 1190 do ip = 1, npoints 1191 in_box(ip) = sum(xx(1:sb%dim, ip)**2) <= (sb%rsize+DELTA)**2 1192 end do 1193 1194 case(CYLINDER) 1195 do ip = 1, npoints 1196 rr = sqrt(sum(xx(2:sb%dim, ip)**2)) 1197 in_box(ip) = rr <= sb%rsize + DELTA 1198 if(sb%periodic_dim >= 1) then 1199 in_box(ip) = in_box(ip) .and. xx(1, ip) >= -sb%xsize - DELTA 1200 in_box(ip) = in_box(ip) .and. xx(1, ip) <= sb%xsize - DELTA 1201 else 1202 in_box(ip) = in_box(ip) .and. abs(xx(1, ip)) <= sb%xsize + DELTA 1203 end if 1204 end do 1205 1206 case(MINIMUM) 1207 1208 if(sb%rsize > M_ZERO) then 1209 radius = sb%rsize 1210 else 1211 radius = M_ZERO 1212 do iatom = 1, geo%natoms 1213 if(species_def_rsize(geo%atom(iatom)%species) < -M_EPSILON) then 1214 write(message(1),'(a,a,a)') 'Using default radii for minimum box, but radius for ', & 1215 trim(species_label(geo%atom(iatom)%species)), ' is negative or undefined.' 1216 message(2) = "Define it properly in the Species block or set the Radius variable explicitly." 1217 call messages_fatal(2, namespace=namespace) 1218 end if 1219 radius = max(radius, species_def_rsize(geo%atom(iatom)%species)) 1220 end do 1221 end if 1222 1223 radius = radius + DELTA 1224 1225 SAFE_ALLOCATE(nlist(1:npoints)) 1226 1227 if(sb%rsize > M_ZERO) then 1228 nullify(list) 1229 call lookup_get_list(sb%atom_lookup, npoints, xx, radius, nlist) 1230 else 1231 call lookup_get_list(sb%atom_lookup, npoints, xx, radius, nlist, list = list) 1232 end if 1233 1234 if(sb%rsize > M_ZERO) then 1235 do ip = 1, npoints 1236 in_box(ip) = (nlist(ip) /= 0) 1237 end do 1238 else 1239 do ip = 1, npoints 1240 in_box(ip) = .false. 1241 do ilist = 1, nlist(ip) 1242 iatom = list(ilist, ip) 1243 dist2 = sum((xx(1:sb%dim, ip) - geo%atom(iatom)%x(1:sb%dim))**2) 1244 if(dist2 < species_def_rsize(geo%atom(iatom)%species)**2) then 1245 in_box(ip) = .true. 1246 exit 1247 end if 1248 end do 1249 end do 1250 end if 1251 1252 SAFE_DEALLOCATE_A(nlist) 1253 SAFE_DEALLOCATE_P(list) 1254 1255 case(PARALLELEPIPED, HYPERCUBE) 1256 llimit(1:sb%dim) = -sb%lsize(1:sb%dim) - DELTA 1257 ulimit(1:sb%dim) = sb%lsize(1:sb%dim) + DELTA 1258 ulimit(1:sb%periodic_dim) = sb%lsize(1:sb%periodic_dim) - DELTA 1259 1260 do ip = 1, npoints 1261 in_box(ip) = all(xx(1:sb%dim, ip) >= llimit(1:sb%dim) .and. xx(1:sb%dim, ip) <= ulimit(1:sb%dim)) 1262 end do 1263 1264#if defined(HAVE_GDLIB) 1265! Why the minus sign for y? Explanation: http://biolinx.bios.niu.edu/bios546/gd_mod.htm 1266! For reasons that probably made sense to someone at some time, computer graphic coordinates are not the same 1267! as in standard graphing. ... The top left corner of the screen is (0,0). 1268 1269 case(BOX_IMAGE) 1270 do ip = 1, npoints 1271 ix = nint(( xx(1, ip) + sb%lsize(1)) * sb%image_size(1) / (M_TWO * sb%lsize(1))) 1272 iy = nint((-xx(2, ip) + sb%lsize(2)) * sb%image_size(2) / (M_TWO * sb%lsize(2))) 1273 call gdlib_image_get_pixel_rgb(sb%image, ix, iy, red, green, blue) 1274 in_box(ip) = (red == 255) .and. (green == 255) .and. (blue == 255) 1275 end do 1276#endif 1277 1278 case(BOX_USDEF) 1279 ! is it inside the user-given boundaries? 1280 do ip = 1, npoints 1281 in_box(ip) = all(xx(1:sb%dim, ip) >= -sb%lsize(1:sb%dim) - DELTA) & 1282 .and. all(xx(1:sb%dim, ip) <= sb%lsize(1:sb%dim) + DELTA) 1283 1284 ! and inside the simulation box? 1285 do idir = 1, sb%dim 1286 xx(idir, ip) = units_from_atomic(units_inp%length, xx(idir, ip)) 1287 end do 1288 rr = sqrt(sum(xx(1:sb%dim, ip)**2)) 1289 call parse_expression(re, im, sb%dim, xx(:, ip), rr, M_ZERO, sb%user_def) 1290 in_box(ip) = in_box(ip) .and. (re /= M_ZERO) 1291 end do 1292 end select 1293 1294 SAFE_DEALLOCATE_A(xx) 1295 1296 end subroutine simul_box_in_box_vec 1297 1298 1299 !-------------------------------------------------------------- 1300 logical pure function simul_box_is_periodic(sb) 1301 type(simul_box_t), intent(in) :: sb 1302 1303 simul_box_is_periodic = sb%periodic_dim > 0 1304 1305 end function simul_box_is_periodic 1306 1307 1308 !-------------------------------------------------------------- 1309 logical pure function simul_box_has_zero_bc(sb) 1310 type(simul_box_t), intent(in) :: sb 1311 1312 simul_box_has_zero_bc = .not. simul_box_is_periodic(sb) 1313 1314 end function simul_box_has_zero_bc 1315 1316 1317 !-------------------------------------------------------------- 1318 subroutine simul_box_dump(sb, namespace, dir, filename, mpi_grp, ierr) 1319 type(simul_box_t), intent(in) :: sb 1320 type(namespace_t), intent(in) :: namespace 1321 character(len=*), intent(in) :: dir 1322 character(len=*), intent(in) :: filename 1323 type(mpi_grp_t), intent(in) :: mpi_grp 1324 integer, intent(out) :: ierr 1325 1326 integer :: iunit, idir 1327 1328 PUSH_SUB(simul_box_dump) 1329 1330 ierr = 0 1331 1332 iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='write', & 1333 position="append", die=.false., grp=mpi_grp) 1334 if (iunit <= 0) then 1335 ierr = ierr + 1 1336 message(1) = "Unable to open file '"//trim(dir)//"/"//trim(filename)//"'." 1337 call messages_warning(1, namespace=namespace) 1338 else 1339 !Only root writes 1340 if (mpi_grp_is_root(mpi_grp)) then 1341 write(iunit, '(a)') dump_tag 1342 write(iunit, '(a20,i4)') 'box_shape= ', sb%box_shape 1343 write(iunit, '(a20,i4)') 'dim= ', sb%dim 1344 write(iunit, '(a20,i4)') 'periodic_dim= ', sb%periodic_dim 1345 write(iunit, '(a20,i4)') 'transport_dim= ', 0 ! sb%transport_dim 1346 select case(sb%box_shape) 1347 case(SPHERE, MINIMUM) 1348 write(iunit, '(a20,e22.14)') 'rsize= ', sb%rsize 1349 write(iunit, '(a20,99e22.14)') 'lsize= ', sb%lsize(1:sb%dim) 1350 case(CYLINDER) 1351 write(iunit, '(a20,e22.14)') 'rsize= ', sb%rsize 1352 write(iunit, '(a20,e22.14)') 'xlength= ', sb%xsize 1353 write(iunit, '(a20,99e22.14)') 'lsize= ', sb%lsize(1:sb%dim) 1354 case(PARALLELEPIPED) 1355 write(iunit, '(a20,99e22.14)') 'lsize= ', sb%lsize(1:sb%dim) 1356 case(BOX_USDEF) 1357 write(iunit, '(a20,99e22.14)') 'lsize= ', sb%lsize(1:sb%dim) 1358 write(iunit, '(a20,a1024)') 'user_def= ', sb%user_def 1359 end select 1360 write(iunit, '(a20,99e22.14)') 'box_offset= ', (M_ZERO, idir = 1, sb%dim) 1361 write(iunit, '(a20,l7)') 'mr_flag= ', sb%mr_flag 1362 if(sb%mr_flag) then 1363 write(iunit, '(a20,i4)') 'num_areas= ',sb%hr_area%num_areas 1364 write(iunit, '(a20,i4)') 'num_radii= ',sb%hr_area%num_radii 1365 do idir = 1, sb%hr_area%num_radii 1366 write(iunit, '(a10,i2.2,a9,e22.14)') 'mr_radius_', idir, '= ',sb%hr_area%radius(idir) 1367 end do 1368 do idir = 1, sb%dim 1369 write(iunit, '(a7,i1,a13,e22.14)') 'center(', idir, ')= ',sb%hr_area%center(idir) 1370 end do 1371 end if 1372 do idir = 1, sb%dim 1373 write(iunit, '(a9,i1,a11,99e22.14)') 'rlattice(', idir, ')= ', & 1374 sb%rlattice_primitive(1:sb%dim, idir) 1375 end do 1376 end if 1377 1378 call io_close(iunit, grp=mpi_grp) 1379 end if 1380 1381 POP_SUB(simul_box_dump) 1382 end subroutine simul_box_dump 1383 1384 1385 ! -------------------------------------------------------------- 1386 subroutine simul_box_load(sb, namespace, dir, filename, mpi_grp, ierr) 1387 type(simul_box_t), intent(inout) :: sb 1388 type(namespace_t), intent(in) :: namespace 1389 character(len=*), intent(in) :: dir 1390 character(len=*), intent(in) :: filename 1391 type(mpi_grp_t), intent(in) :: mpi_grp 1392 integer, intent(out) :: ierr 1393 1394 integer :: iunit, idim, il, err 1395 character(len=20) :: str 1396 character(len=100), allocatable :: lines(:) 1397 FLOAT :: rlattice_primitive(1:MAX_DIM, 1:MAX_DIM) 1398 1399 PUSH_SUB(simul_box_load) 1400 1401 ierr = 0 1402 1403 iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='read', & 1404 status="old", die=.false., grp=mpi_grp) 1405 if (iunit <= 0) then 1406 ierr = ierr + 1 1407 message(1) = "Unable to open file '"//trim(dir)//"/"//trim(filename)//"'." 1408 call messages_warning(1, namespace=namespace) 1409 else 1410 ! Find the dump tag. 1411 call iopar_find_line(mpi_grp, iunit, dump_tag, err) 1412 if (err /= 0) ierr = ierr + 2 1413 1414 SAFE_ALLOCATE(lines(1:4)) 1415 if (ierr == 0) then 1416 call iopar_read(mpi_grp, iunit, lines, 4, err) 1417 if (err == 0) then 1418 read(lines(1), *) str, sb%box_shape 1419 read(lines(2), *) str, sb%dim 1420 read(lines(3), *) str, sb%periodic_dim 1421 read(lines(4), *) str, il ! sb%transport_dim 1422 else 1423 ierr = ierr + 2**2 1424 end if 1425 1426 select case (sb%box_shape) 1427 case(SPHERE, MINIMUM) 1428 call iopar_read(mpi_grp, iunit, lines, 2, err) 1429 if (err /= 0) then 1430 ierr = ierr + 2**3 1431 else 1432 read(lines(1), *) str, sb%rsize 1433 read(lines(2), *) str, sb%lsize(1:sb%dim) 1434 end if 1435 case(CYLINDER) 1436 call iopar_read(mpi_grp, iunit, lines, 3, err) 1437 if (err /= 0) then 1438 ierr = ierr + 2**4 1439 else 1440 read(lines(1), *) str, sb%rsize 1441 read(lines(2), *) str, sb%xsize 1442 read(lines(3), *) str, sb%lsize(1:sb%dim) 1443 end if 1444 case(PARALLELEPIPED) 1445 call iopar_read(mpi_grp, iunit, lines, 1, err) 1446 if (err /= 0) then 1447 ierr = ierr + 2**5 1448 else 1449 read(lines(1), *) str, sb%lsize(1:sb%dim) 1450 end if 1451 case(BOX_USDEF) 1452 call iopar_read(mpi_grp, iunit, lines, 2, err) 1453 if (err /= 0) then 1454 ierr = ierr + 2**6 1455 else 1456 read(lines(1), *) str, sb%lsize(1:sb%dim) 1457 read(lines(2), *) str, sb%user_def 1458 end if 1459 end select 1460 1461 sb%mr_flag = .false. 1462 call iopar_read(mpi_grp, iunit, lines, 2, err) 1463 if (err /= 0) then 1464 ierr = ierr + 2**7 1465 else 1466 ! lines(1) was sb%box_offset, now removed 1467 read(lines(2),'(a20,l7)') str, sb%mr_flag 1468 end if 1469 1470 SAFE_DEALLOCATE_A(lines) 1471 1472 if (sb%mr_flag) then 1473 SAFE_ALLOCATE(lines(1:2)) 1474 call iopar_read(mpi_grp, iunit, lines, 2, err) 1475 if (err /= 0) then 1476 ierr = ierr + 2**8 1477 else 1478 read(lines(1),*) str, sb%hr_area%num_areas 1479 read(lines(2),*) str, sb%hr_area%num_radii 1480 end if 1481 SAFE_DEALLOCATE_A(lines) 1482 1483 SAFE_ALLOCATE(sb%hr_area%radius(1:sb%hr_area%num_radii)) 1484 SAFE_ALLOCATE(lines(1:sb%hr_area%num_radii)) 1485 sb%hr_area%num_radii = 0 1486 call iopar_read(mpi_grp, iunit, lines, sb%hr_area%num_radii, err) 1487 if (err /= 0) then 1488 ierr = ierr + 2**9 1489 else 1490 do il = 1, sb%hr_area%num_radii 1491 read(lines(1),*) str, sb%hr_area%radius(il) 1492 end do 1493 end if 1494 SAFE_DEALLOCATE_A(lines) 1495 1496 SAFE_ALLOCATE(lines(1:sb%dim)) 1497 call iopar_read(mpi_grp, iunit, lines, sb%dim, err) 1498 if (err /= 0) then 1499 ierr = ierr + 2**10 1500 else 1501 do idim = 1, sb%dim 1502 read(lines(1), *) str, sb%hr_area%center(idim) 1503 end do 1504 end if 1505 SAFE_DEALLOCATE_A(lines) 1506 end if 1507 1508 1509 SAFE_ALLOCATE(lines(1:sb%dim)) 1510 call iopar_read(mpi_grp, iunit, lines, sb%dim, err) 1511 if (err /= 0) then 1512 ierr = ierr + 2**11 1513 else 1514 do idim = 1, sb%dim 1515 read(lines(idim), *) str, rlattice_primitive(1:sb%dim, idim) 1516 end do 1517 end if 1518 SAFE_DEALLOCATE_A(lines) 1519 1520 end if 1521 1522 call io_close(iunit, grp=mpi_grp) 1523 end if 1524 1525 if (ierr == 0) then 1526 call simul_box_build_lattice(sb, namespace, rlattice_primitive) 1527 end if 1528 1529 POP_SUB(simul_box_load) 1530 end subroutine simul_box_load 1531 1532 1533 ! -------------------------------------------------------------- 1534 recursive subroutine simul_box_copy(sbout, sbin) 1535 type(simul_box_t), intent(out) :: sbout 1536 type(simul_box_t), intent(in) :: sbin 1537 1538 PUSH_SUB(simul_box_copy) 1539 1540 sbout%box_shape = sbin%box_shape 1541 sbout%rsize = sbin%rsize 1542 sbout%xsize = sbin%xsize 1543 sbout%lsize = sbin%lsize 1544 sbout%image = sbin%image 1545 sbout%user_def = sbin%user_def 1546 sbout%rlattice = sbin%rlattice 1547 sbout%rlattice_primitive = sbin%rlattice_primitive 1548 sbout%klattice = sbin%klattice 1549 sbout%klattice_primitive = sbin%klattice_primitive 1550 sbout%volume_element = sbin%volume_element 1551 sbout%dim = sbin%dim 1552 sbout%periodic_dim = sbin%periodic_dim 1553 sbout%mr_flag = sbin%mr_flag 1554 sbout%hr_area%num_areas = sbin%hr_area%num_areas 1555 sbout%hr_area%num_radii = sbin%hr_area%num_radii 1556 sbout%hr_area%center(1:sbin%dim)=sbin%hr_area%center(1:sbin%dim) 1557 1558 call kpoints_copy(sbin%kpoints, sbout%kpoints) 1559 1560 if(sbout%mr_flag) then 1561 SAFE_ALLOCATE(sbout%hr_area%radius(1:sbout%hr_area%num_radii)) 1562 sbout%hr_area%radius(1:sbout%hr_area%num_radii) = sbin%hr_area%radius(1:sbout%hr_area%num_radii) 1563 end if 1564 1565 call lookup_copy(sbin%atom_lookup, sbout%atom_lookup) 1566 1567 if(simul_box_is_periodic(sbin)) call symmetries_copy(sbin%symm, sbout%symm) 1568 1569 POP_SUB(simul_box_copy) 1570 end subroutine simul_box_copy 1571 1572 ! ----------------------------------------------------- 1573 1574 subroutine simul_box_check_atoms_are_too_close(geo, sb, namespace) 1575 type(geometry_t), intent(in) :: geo 1576 type(simul_box_t), intent(in) :: sb 1577 type(namespace_t), intent(in) :: namespace 1578 1579 FLOAT :: mindist 1580 FLOAT, parameter :: threshold = CNST(1e-5) 1581 1582 PUSH_SUB(simul_box_check_atoms_are_too_close) 1583 1584 if(geo%natoms == 1) then 1585 POP_SUB(simul_box_check_atoms_are_too_close) 1586 return 1587 end if 1588 1589 mindist = simul_box_min_distance(geo, sb, real_atoms_only = .false.) 1590 if(mindist < threshold) then 1591 write(message(1), '(a)') "Some of the atoms seem to sit too close to each other." 1592 write(message(2), '(a)') "Please review your input files and the output geometry (in 'static/')." 1593 write(message(3), '(a, f12.6, 1x, a)') "Minimum distance = ", & 1594 units_from_atomic(units_out%length, mindist), trim(units_abbrev(units_out%length)) 1595 call messages_warning(3, namespace=namespace) 1596 1597 ! then write out the geometry, whether asked for or not in Output variable 1598 call io_mkdir(STATIC_DIR, namespace) 1599 call geometry_write_xyz(geo, trim(STATIC_DIR)//'/geometry', namespace) 1600 end if 1601 1602 if(simul_box_min_distance(geo, sb, real_atoms_only = .true.) < threshold) then 1603 message(1) = "It cannot be correct to run with physical atoms so close." 1604 call messages_fatal(1, namespace=namespace) 1605 end if 1606 1607 POP_SUB(simul_box_check_atoms_are_too_close) 1608 end subroutine simul_box_check_atoms_are_too_close 1609 1610 ! --------------------------------------------------------- 1611 FLOAT function simul_box_min_distance(geo, sb, real_atoms_only) result(rmin) 1612 type(geometry_t), intent(in) :: geo 1613 type(simul_box_t), intent(in) :: sb 1614 logical, optional, intent(in) :: real_atoms_only 1615 1616 integer :: iatom, jatom, idir 1617 FLOAT :: xx(MAX_DIM) 1618 logical :: real_atoms_only_ 1619 type(species_t), pointer :: species 1620 1621 PUSH_SUB(simul_box_min_distance) 1622 1623 real_atoms_only_ = optional_default(real_atoms_only, .false.) 1624 1625 rmin = huge(rmin) 1626 do iatom = 1, geo%natoms 1627 call atom_get_species(geo%atom(iatom), species) 1628 if(real_atoms_only_ .and. .not. species_represents_real_atom(species)) cycle 1629 do jatom = iatom + 1, geo%natoms 1630 call atom_get_species(geo%atom(iatom), species) 1631 if(real_atoms_only_ .and. .not. species_represents_real_atom(species)) cycle 1632 xx(:) = abs(geo%atom(iatom)%x(:) - geo%atom(jatom)%x(:)) 1633 do idir = 1, sb%periodic_dim 1634 xx(idir) = xx(idir) - M_TWO * sb%lsize(idir) * floor(xx(idir)/(M_TWO * sb%lsize(idir)) + M_HALF) 1635 end do 1636 rmin = min(sqrt(sum(xx**2)), rmin) 1637 end do 1638 end do 1639 1640 if(.not. (geo%only_user_def .and. real_atoms_only_)) then 1641 ! what if the nearest neighbors are periodic images? 1642 do idir = 1, sb%periodic_dim 1643 rmin = min(rmin, abs(sb%lsize(idir))) 1644 end do 1645 end if 1646 1647 POP_SUB(simul_box_min_distance) 1648 end function simul_box_min_distance 1649 1650 1651 ! --------------------------------------------------------- 1652 subroutine simul_box_symmetry_check(this, geo, dim, namespace) 1653 type(simul_box_t), intent(in) :: this 1654 type(geometry_t), intent(in) :: geo 1655 integer, intent(in) :: dim 1656 type(namespace_t), intent(in) :: namespace 1657 1658 integer :: iop, iatom, iatom_symm 1659 FLOAT :: ratom(1:MAX_DIM) 1660 1661 PUSH_SUB(simul_box_symmetry_check) 1662 1663 ! We want to use for instance that 1664 ! 1665 ! \int dr f(Rr) V_iatom(r) \nabla f(R(v)) = R\int dr f(r) V_iatom(R*r) f(r) 1666 ! 1667 ! and that the operator R should map the position of atom 1668 ! iatom to the position of some other atom iatom_symm, so that 1669 ! 1670 ! V_iatom(R*r) = V_iatom_symm(r) 1671 ! 1672 do iop = 1, symmetries_number(this%symm) 1673 if(iop == symmetries_identity_index(this%symm)) cycle 1674 1675 do iatom = 1, geo%natoms 1676 ratom = M_ZERO 1677 if(geo%reduced_coordinates) then 1678 ratom(1:this%dim) = symm_op_apply_red(this%symm%ops(iop), geo%atom(iatom)%x) 1679 else 1680 ratom(1:this%dim) = symm_op_apply_cart(this%symm%ops(iop), geo%atom(iatom)%x) 1681 end if 1682 1683 call simul_box_periodic_atom_in_box(this, geo, ratom) 1684 1685 ! find iatom_symm 1686 do iatom_symm = 1, geo%natoms 1687 if(all(abs(ratom(1:dim) - geo%atom(iatom_symm)%x(1:dim)) < CNST(1.0e-5))) exit 1688 end do 1689 1690 if(iatom_symm > geo%natoms) then 1691 write(message(1),'(a,i6)') 'Internal error: could not find symmetric partner for atom number', iatom 1692 write(message(2),'(a,i3,a)') 'with symmetry operation number ', iop, '.' 1693 call messages_fatal(2, namespace=namespace) 1694 end if 1695 1696 end do 1697 end do 1698 1699 POP_SUB(simul_box_symmetry_check) 1700 end subroutine simul_box_symmetry_check 1701 1702end module simul_box_oct_m 1703 1704!! Local Variables: 1705!! mode: f90 1706!! coding: utf-8 1707!! End: 1708