1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 9! This code has been re-coded by: 10! Nick Papior Andersen to handle other than standard constrainst. 11 12! Using CG/Broyden optimization *should* work as it is done per 13! force magnitude, and we correct for that. 14 15module m_fixed 16 17 use precision 18 19 implicit none 20 21 public :: init_fixed, resetFixedPointers 22 public :: print_fixed 23 public :: fixed 24 25 ! Allows to check whether certain 26 ! atoms are fixed 27 public :: is_fixed, is_constr 28 29 private 30 31 integer, parameter :: TYPE_LEN = 20 32 type tFix 33 ! Number of atoms belonging to this fix 34 integer :: n = 0 35 ! the type of fixation 36 character(TYPE_LEN) :: type 37 ! Atoms in this fixation 38 integer, allocatable :: a(:) 39 ! direction of fixation 40 ! Note that this is a single value for all 41 ! If more atoms belong to one fixation 42 ! we fix them all to move along that direction 43 real(dp) :: fix(3) = 0._dp 44 end type tFix 45 46 ! Containers for fixations 47 integer, save :: N_fix = 0 48 type(tFix), allocatable, save :: fixs(:) 49 50 ! Stress fixation 51 real(dp), save :: xs(6) = 0._dp 52 53 ! Use constr routine 54 logical, save :: use_constr = .false. 55 56 ! Fix angles 57 logical, save :: cell_angle(3) = .false. ! alpha, beta, gamma 58 ! Fix cell-vectors 59 logical, save :: cell_vector(3) = .false. ! A1, A2, A3 60 61contains 62 63 subroutine resetFixedPointers( ) 64 65 integer :: i 66 67 if ( N_fix == 0 ) return 68 69 do i = 1 , N_fix 70 deallocate(fixs(i)%a) 71 enddo 72 deallocate(fixs) 73 N_fix = 0 74 75 end subroutine resetFixedPointers 76 77 ! ********************************************************************** 78 ! Reads and imposes required constraints to atomic displacements by 79 ! making zero the forces in those directions. Constraints are specified 80 ! by the FDF data block GeometryConstraints (see example below). 81 ! Only types position and routine implemented in this version. 82 ! Written by J.M.Soler. Feb., 1998 83 ! Modified by P. Ordejon to output the total number of constraints 84 ! imposed. June, 2003 85 ! Rewritten by Nick Papior for additional constraint types and 86 ! printing of explicit constraints. 87 ! January, 2015 88 ! *********** INPUT **************************************************** 89 ! real*8 cell(3,3) : Lattice vectors 90 ! real*8 stress( 3,3) : Stress tensor 91 ! integer na : Number of atoms 92 ! integer isa(na) : Species indexes 93 ! real*8 amass(na) : Atomic masses 94 ! real*8 xa(3,na) : Atomic cartesian coordinates 95 ! real*8 fa(3,na) : Atomic forces 96 ! *********** OUTPUT *************************************************** 97 ! real*8 cstress( 3,3) : Constrained stress tensor 98 ! real*8 cfa(3,na) : Constrained atomic forces 99 ! integer ntcon : Total number of position constraints imposed 100 ! *********** UNITS **************************************************** 101 ! Units are arbitrary but cell and xa must be in the same units 102 ! *********** BEHAVIOUR ************************************************ 103 ! cstress may be the same physical array as stress, and cfa the same 104 ! as fa, i.e. it is allowed: 105 ! call fixed( cell, stress, na, isa, amass, xa, fa, stress, fa, ntcon ) 106 ! NOTE: This is NOT good practice, and might be flagged by some compilers 107 ! *********** USAGE **************************************************** 108 ! Example: consider a diatomic molecule (atoms 1 and 2) above a surface, 109 ! represented by a slab of 5 atomic layers, with 10 atoms per layer. 110 ! To fix the cell height, the slab's botom layer (last 10 atoms), 111 ! the molecule's interatomic distance, its height above the surface 112 ! (but not its azimutal orientation and lateral position), and the 113 ! relative height of the two atoms: 114 ! 115 ! %block GeometryConstraints 116 ! cellside c 117 ! cellangle alpha beta gamma 118 ! position from -1 to -10 119 ! rigid 1 2 120 ! center 1 2 0.0 0.0 1.0 121 ! routine constr 122 ! stress 1 2 3 123 ! %endblock GeometryConstraints 124 ! 125 ! where constr is the following user-written subroutine: 126 ! 127 ! subroutine constr( cell, na, isa, amass, xa, stress, fa, ntcon ) 128 !c real*8 cell(3,3) : input lattice vectors (Bohr) 129 !c integer na : input number of atoms 130 !c integer isa(na) : input species indexes 131 !c real*8 amass(na) : input atomic masses 132 !c real*8 xa(3,na) : input atomic cartesian coordinates (Bohr) 133 !c real*8 stress( 3,3) : input/output stress tensor (Ry/Bohr**3) 134 !c real*8 fa(3,na) : input/output atomic forces (Ry/Bohr) 135 !c integer ntcon : total number of position constraints, accumulative 136 ! integer na, isa(na), ntcon 137 ! double precision amass(na), cell(3,3), fa(3,na), 138 ! . stress(3,3), xa(3,na), fz 139 ! fz = fa(3,1) + fa(3,2) 140 ! fa(3,1) = fz*amass(1)/(amass(1)+amass(2)) 141 ! fa(3,2) = fz*amass(2)/(amass(1)+amass(2)) 142 ! ntcon = ntcon+1 143 ! end 144 ! ********************************************************************** 145 146 subroutine fixed( cell, stress, na, isa, amass, xa, fa, cstress, cfa, ntcon , & 147 magnitude_usage ) 148 149 use intrinsic_missing, only : VNORM, VEC_PROJ 150 151 integer, intent(in) :: na, isa(na) 152 real(dp), intent(in) :: amass(na), cell(3,3), fa(3,na), stress(3,3), xa(3,na) 153 integer, intent(out) :: ntcon 154 real(dp), intent(out) :: cfa(3,na), cstress(3,3) 155 156 ! Whether we use a routine to find the gradient 157 ! of the magnitudes, in which case we do not 158 ! want to scale with mass for certain optimizations. 159 logical, intent(in), optional :: magnitude_usage 160 161 integer :: ix, jx, ia 162 integer :: if, N, i 163 character(len=TYPE_LEN) :: namec 164 real(dp) :: ca(3), am, cf(3) 165 logical :: lmag_use 166 167#ifdef DEBUG 168 call write_debug( ' PRE fixed' ) 169#endif 170 171 ! Copy stress and forces to output arrays 172 do ix = 1,3 173 do jx = 1,3 174 cstress(jx,ix) = stress(jx,ix) 175 end do 176 end do 177 do ia = 1,na 178 do ix = 1,3 179 cfa(ix,ia) = fa(ix,ia) 180 end do 181 end do 182 183 ! Initialize the number of relaxed degrees of freedom 184 ntcon = 0 185 186 ! If we should call the routine 187 if ( use_constr ) then 188 call constr( cell, na, isa, amass, xa, cstress, cfa, ntcon ) 189 end if 190 191 ! apply the stress constraint 192 cstress(1,1) = cstress(1,1) - xs(1) * cstress(1,1) 193 cstress(2,2) = cstress(2,2) - xs(2) * cstress(2,2) 194 cstress(3,3) = cstress(3,3) - xs(3) * cstress(3,3) 195 cstress(3,2) = cstress(3,2) - xs(4) * cstress(3,2) 196 cstress(2,3) = cstress(2,3) - xs(4) * cstress(2,3) 197 cstress(3,1) = cstress(3,1) - xs(5) * cstress(3,1) 198 cstress(1,3) = cstress(1,3) - xs(5) * cstress(1,3) 199 cstress(1,2) = cstress(1,2) - xs(6) * cstress(1,2) 200 cstress(2,1) = cstress(2,1) - xs(6) * cstress(2,1) 201 202 ! stress-directions will be symmetrized 203 ! This constrain is equivalent to only 204 ! allow the same expansion factor in front 205 ! of each cell-vector. 206 if ( cell_angle(1) ) then 207 ! alpha angle, ang(B,C) 208 am = ( cstress(2,2) + cstress(3,3) ) * .5_dp 209 ! Kill everything else... 210 cstress(:,:) = 0._dp 211 cstress(2,2) = am 212 cstress(3,3) = am 213 end if 214 if ( cell_angle(2) ) then 215 ! beta angle, ang(A,C) 216 am = ( cstress(1,1) + cstress(3,3) ) * .5_dp 217 cstress(:,:) = 0._dp 218 cstress(1,1) = am 219 cstress(3,3) = am 220 end if 221 if ( cell_angle(3) ) then 222 ! gamma angle, ang(A,B) 223 am = ( cstress(1,1) + cstress(2,2) ) * .5_dp 224 cstress(:,:) = 0._dp 225 cstress(1,1) = am 226 cstress(2,2) = am 227 end if 228 229 ! should be done after angles... 230 do if = 1 , 3 231 if ( cell_vector(if) ) then 232 cstress(:,if) = 0._dp 233 cstress(if,:) = 0._dp 234 end if 235 end do 236 237 lmag_use = .false. 238 if ( present(magnitude_usage) ) lmag_use = magnitude_usage 239 240 241 if ( N_fix > 0 ) then 242 do if = 1 , N_fix 243 244 ! apply the designated constraints 245 246 N = fixs(if)%N 247 248 namec = fixs(if)%type 249 250 ! Select type of constraint 251 if ( namec == 'pos' ) then 252 253 ! this is the easy one, all atoms should not move 254 do i = 1 , N 255 256 ia = fixs(if)%a(i) 257 cfa(:,ia) = 0._dp 258 259 end do 260 261 ! we remove N * 3 degrees of freedom 262 ntcon = ntcon + 3 * N 263 264 else if ( namec == 'pos-dir' ) then 265 266 do i = 1 , N 267 268 ia = fixs(if)%a(i) 269 270 ! Calculate the force projected onto the fixation vector 271 cfa(:,ia) = cfa(:,ia) - VEC_PROJ( fixs(if)%fix , cfa(:,ia) ) 272 273 end do 274 275 ! Only one direction is fixed, hence all atoms 276 ! remove one degree of freedom 277 ntcon = ntcon + N 278 279 else if ( namec == 'center' ) then 280 281 ! Maintain the same center of the molecule 282 cf(:) = 0._dp 283 284 ! we center the molecule by constraining the 285 ! average acceleration to 0 (for MD) 286 ! or by subtracting the average force for magnitude used forces 287 do i = 1 , N 288 289 ia = fixs(if)%a(i) 290 291 if ( lmag_use ) then 292 cf(:) = cf(:) + cfa(:,ia) 293 else 294 cf(:) = cf(:) + cfa(:,ia) / amass(ia) 295 end if 296 297 end do 298 299 ! This is the average force/acceleration 300 cf(:) = cf(:) / real(N,dp) 301 302 ! fix for the correct direction, project onto the fixation vector 303 do i = 1 , N 304 305 ia = fixs(if)%a(i) 306 307 if ( lmag_use ) then 308 ! subtract the average force so that the net-force is zero. 309 cfa(:,ia) = cfa(:,ia) - cf(:) 310 else 311 cfa(:,ia) = cfa(:,ia) - cf(:) * amass(ia) 312 end if 313 314 end do 315 316 ! We constrain the translational part which means that 317 ! we constrain 3 degrees of freedom 318 ntcon = ntcon + 3 319 320 else if ( namec == 'center-dir' ) then 321 322 ! Maintain the same center of the molecule 323 cf(:) = 0._dp 324 325 ! we center the molecule by constraining the 326 ! average acceleration to 0 (for MD) 327 ! or by subtracting the average force for magnitude used forces 328 do i = 1 , N 329 330 ia = fixs(if)%a(i) 331 332 if ( lmag_use ) then 333 cf(:) = cf(:) + cfa(:,ia) 334 else 335 cf(:) = cf(:) + cfa(:,ia) / amass(ia) 336 end if 337 338 end do 339 340 ! This is the average force/acceleration 341 cf(:) = cf(:) / real(N,dp) 342 cf(:) = cf(:) - VEC_PROJ( fixs(if)%fix, cf ) 343 344 ! fix for the correct direction, project onto the fixation vector 345 do i = 1 , N 346 347 ia = fixs(if)%a(i) 348 349 if ( lmag_use ) then 350 ! subtract the average force so that the net-force is zero. 351 cfa(:,ia) = cfa(:,ia) - cf(:) 352 else 353 cfa(:,ia) = cfa(:,ia) - cf(:) * amass(ia) 354 end if 355 356 end do 357 358 ! We constrain the translational part along one 359 ! direction, thus we constrain 1 degree of freedom 360 ntcon = ntcon + 1 361 362 else if ( namec == 'com' ) then 363 364 ! Calculate center of mass of the molecule 365 ca(:) = 0._dp 366 am = 0._dp 367 368 do i = 1 , N 369 370 ia = fixs(if)%a(i) 371 372 ! center of mass 373 ca(:) = ca(:) + xa(:,ia) * amass(ia) 374 am = am + amass(ia) 375 376 end do 377 378 ! get the center of mass 379 ca(:) = ca(:) / am 380 381 ! correct the forces so that the center of mass is the same 382 383 call die('Center of mass does not currently work') 384 385 else if ( namec == 'rigid' ) then 386 387 ! Calculate total force on the molecule 388 ! this is done using the center-of-force method 389 cf(:) = 0._dp 390 am = 0._dp 391 392 do i = 1 , N 393 394 ia = fixs(if)%a(i) 395 396 ! sum the force 397 cf(:) = cf(:) + cfa(:,ia) 398 am = am + amass(ia) 399 400 end do 401 402 if ( lmag_use ) then 403 ! use the average mass. 404 cf(:) = cf(:) / real(N,dp) 405 else 406 cf(:) = cf(:) / am 407 end if 408 409 do i = 1 , N 410 411 ia = fixs(if)%a(i) 412 413 if ( lmag_use ) then 414 cfa(:,ia) = cf(:) 415 else 416 cfa(:,ia) = cf(:) * amass(ia) 417 end if 418 419 end do 420 421 ! we constrain (N - 1) * 3 degrees of freedom 422 ! (only the relative positions are constrained) 423 ntcon = ntcon + ( N - 1 ) * 3 424 425 else if ( namec == 'rigid-dir' ) then 426 427 ! Calculate total force on the molecule 428 ! this is done using the center-of-force method 429 cf(:) = 0._dp 430 am = 0._dp 431 432 ! not so straight forward constraint 433 do i = 1 , N 434 435 ia = fixs(if)%a(i) 436 437 ! sum the force 438 cf(:) = cf(:) + cfa(:,ia) 439 am = am + amass(ia) 440 441 end do 442 443 if ( lmag_use ) then 444 ! use the average force 445 cf(:) = cf(:) / real(N,dp) 446 else 447 ! use the average mass 448 cf(:) = cf(:) / am 449 end if 450 451 ! only set the molecular force in the 452 ! specified direction, project onto the fixation vector 453 cf(:) = cf(:) - VEC_PROJ( fixs(if)%fix(:) , cf(:) ) 454 455 do i = 1 , N 456 457 ia = fixs(if)%a(i) 458 459 if ( lmag_use ) then 460 cfa(:,ia) = cf(:) 461 else 462 cfa(:,ia) = cf(:) * amass(ia) 463 end if 464 465 end do 466 467 ! we constrain N atoms to not move in one direction 468 ! hence we remove N degrees of freedom 469 ntcon = ntcon + N 470 471 else if ( namec == 'rigid-max' ) then 472 473 ! Calculate total force on the molecule 474 ! this is done using the center-of-force method 475 cf(:) = 0._dp 476 am = 0._dp 477 478 do i = 1 , N 479 480 ia = fixs(if)%a(i) 481 482 if ( VNORM(cfa(:,ia)) > VNORM(cf) ) then 483 cf(:) = cfa(:,ia) 484 am = amass(ia) 485 end if 486 487 end do 488 489 if ( lmag_use ) then 490 ! do nothing, we have the max force 491 else 492 cf(:) = cf(:) / am 493 end if 494 495 do i = 1 , N 496 497 ia = fixs(if)%a(i) 498 499 if ( lmag_use ) then 500 cfa(:,ia) = cf(:) 501 else 502 cfa(:,ia) = cf(:) * amass(ia) 503 end if 504 505 end do 506 507 ! we constrain (N - 1) * 3 degrees of freedom (only the relative positions are 508 ! constrained) 509 ntcon = ntcon + ( N - 1 ) * 3 510 511 else if ( namec == 'rigid-max-dir' ) then 512 513 ! Calculate total force on the molecule 514 ! this is done using the center-of-force method 515 cf(:) = 0._dp 516 am = 0._dp 517 518 ! not so straight forward constraint 519 do i = 1 , N 520 521 ia = fixs(if)%a(i) 522 523 if ( VNORM(cfa(:,ia)) > VNORM(cf) ) then 524 cf(:) = cfa(:,ia) 525 am = amass(ia) 526 end if 527 end do 528 529 if ( lmag_use ) then 530 ! do nothing, we have the max force 531 else 532 cf(:) = cf(:) / am 533 end if 534 535 ! only set the molecular force in the 536 ! specified direction, project onto the fixation vector 537 cf(:) = cf(:) - VEC_PROJ( fixs(if)%fix(:) , cf(:) ) 538 539 do i = 1 , N 540 541 ia = fixs(if)%a(i) 542 543 if ( lmag_use ) then 544 cfa(:,ia) = cf(:) 545 else 546 cfa(:,ia) = cf(:) * amass(ia) 547 end if 548 549 end do 550 551 ! we constrain N atoms to not move in one direction 552 ! hence we remove N degrees of freedom 553 ntcon = ntcon + N 554 555 end if 556 557 end do 558 end if 559 560#ifdef DEBUG 561 call write_debug( ' POS fixed' ) 562#endif 563 564 end subroutine fixed 565 566 subroutine print_fixed( ) 567 568 use parallel, only : Node 569 use m_region 570 571 integer :: if 572 character(len=70) :: name 573 574 type(tRgn) :: r 575 576 577 if ( Node /= 0 ) return 578 579 write(*,*) ! newline 580 581 if ( use_constr ) then 582 write(*,'(a)') 'siesta: Constraints using custom constr routine' 583 end if 584 585 if ( any(xs /= 0._dp) ) then 586 write(*,'(a)') 'siesta: Constraint (stress):' 587 write(*,'(3(2(tr2,e11.4),/))') xs(:) 588 end if 589 590 if ( cell_angle(1) ) & 591 write(*,'(a)') 'siesta: Constrain angle alpha (B-C/2-3)' 592 if ( cell_angle(2) ) & 593 write(*,'(a)') 'siesta: Constrain angle beta (A-C/1-3)' 594 if ( cell_angle(3) ) & 595 write(*,'(a)') 'siesta: Constrain angle gamma (A-B/1-2)' 596 597 if ( cell_vector(1) ) & 598 write(*,'(a)') 'siesta: Constraint (vector): A/1' 599 if ( cell_vector(2) ) & 600 write(*,'(a)') 'siesta: Constraint (vector): B/2' 601 if ( cell_vector(3) ) & 602 write(*,'(a)') 'siesta: Constraint (vector): C/3' 603 604 if ( N_fix == 0 ) return 605 606 if ( N_fix > 1 ) then 607 write(*,'(a)') 'siesta: Constraints applied in the following order:' 608 end if 609 610 do if = 1 , N_fix 611 612 call rgn_list(r,fixs(if)%n,fixs(if)%a) 613 r%name = trim(fixs(if)%type) 614 name = 'siesta: Constraint' 615 if ( index(r%name,'dir') > 0 ) then 616 ! the name should include the fixation vector 617 write(name,'(a,2(tr1,f8.5,'',''),tr1,f8.5,a)') & 618 'siesta: Constraint v=[',fixs(if)%fix,']' 619 end if 620 call rgn_print(r, name = trim(name), seq_max = 12) 621 622 end do 623 call rgn_delete(r) 624 625 write(*,*) ! newline 626 627 end subroutine print_fixed 628 629 subroutine init_fixed( cell, na , isa , iza ) 630 631 use fdf 632 use fdf_extra 633 use intrinsic_missing, only : VNORM, VEC_PROJ 634 use parallel, only : IONode 635 636 use m_region 637 638 real(dp), intent(in) :: cell(3,3) 639 integer, intent(in) :: na, isa(na), iza(na) 640 641 ! Internal variables 642 character(len=TYPE_LEN) :: namec 643 644 integer :: sfix, ifix, i, ix, N 645 logical :: add_dir 646 647 type(block_fdf) :: bfdf 648 type(parsed_line), pointer :: pline 649 650 type(tRgn) :: rr, r_tmp 651 652 ! Initialise stress constraints to unconstrained state 653 xs(:) = 0._dp 654 655 ! No fixation atoms 656 N_fix = 0 657 658 ! Look for constraints data block, we also allow another 659 ! constraint block 660 if ( .not. fdf_block('Geometry.Constraints',bfdf) ) return 661 662#ifdef DEBUG 663 call write_debug( ' PRE init_fixed' ) 664#endif 665 666 ! First read in number of constraints 667 do while ( fdf_bline(bfdf,pline) ) 668 669 ! If no names exists, we have no constraint line 670 if ( fdf_bnnames(pline) == 0 ) cycle 671 672 namec = fdf_bnames(pline,1) 673 674 if ( leqi(namec,'position') .or. leqi(namec,'atom') .or. & 675 leqi(namec,'rigid') .or. leqi(namec,'rigid-max') .or. & 676 leqi(namec,'molecule') .or. leqi(namec,'molecule-max') .or. & 677 leqi(namec,'center') .or. & 678 leqi(namec,'center-of-mass') .or. & 679 leqi(namec,'species-i') .or. leqi(namec,'Z') ) then 680 681 ! All these names belong to atomic constraints 682 N_fix = N_fix + 1 683 684 end if 685 686 end do 687 688 ! Allocate space 689 allocate(fixs(N_fix)) 690 691 call fdf_brewind(bfdf) 692 693 ! First read in number of constraints 694 ifix = 0 695 do while ( fdf_bline(bfdf,pline) ) 696 697 ! If no names exists, we have no constraint line 698 if ( fdf_bnnames(pline) == 0 ) cycle 699 700 namec = fdf_bnames(pline,1) 701 ! Indicates whether there should be read in direction 702 ! bound constrictions 703 add_dir = .false. 704 705 ! Take those not specific to atomic positions 706 if ( leqi(namec,'stress') ) then 707 708 N = fdf_bnvalues(pline) 709 710 do i = 1, N 711 712 ix = nint(fdf_bvalues(pline,i)) 713 714 if ( ix < 1 .or. 6 < ix ) then 715 if ( IONOde ) then 716 write(*,'(a)') 'Stress constraint:' 717 write(*,'(6(tr2,i0,'' == '',a,/))') & 718 1,'aa',2,'bb',3,'cc', & 719 4,'bc / cb',5,'ca / ac',6,'ab / ba' 720 end if 721 call die('fixed: Stress restriction not & 722 &with expected input [1:6]') 723 end if 724 725 xs(ix) = 1._dp 726 727 end do 728 729 else if ( leqi(namec,'routine') ) then 730 731 use_constr = .true. 732 733 734 else if ( leqi(namec,'cellangle') .or. leqi(namec,'cell-angle') ) then 735 736 N = fdf_bnnames(pline) 737 738 ! When fixing the cell-angles 739 ! the cell-vectors must be orthogonal 740 ! This is because of non-symmetry in 741 ! the stress-tensor for non-orthogonal 742 ! lattice vectors. 743 744 do i = 2 , N 745 namec = fdf_bnames(pline,i) 746 ! Fix alpha angle 747 if ( leqi(namec,'alpha') ) then 748 if ( is_zero(cell(1,2)) .and. & 749 is_zero(cell(1,3)) ) then 750 cell_angle(1) = .true. 751 else 752 if ( IONOde ) then 753 write(*,'(a)') 'fixed: alpha-angle cannot & 754 &be constrained due to B and C having & 755 &components along Cartesian X-direction.' 756 end if 757 call die('fixed: alpha-angle cannot be constrained.') 758 end if 759 760 else if ( leqi(namec,'beta') ) then 761 if ( is_zero(cell(2,1)) .and. & 762 is_zero(cell(2,3)) ) then 763 cell_angle(2) = .true. 764 else 765 if ( IONOde ) then 766 write(*,'(a)') 'fixed: beta-angle cannot & 767 &be constrained due to A and C having & 768 &components along Cartesian Y-direction.' 769 end if 770 call die('fixed: beta-angle cannot be constrained.') 771 end if 772 773 else if ( leqi(namec,'gamma') ) then 774 if ( is_zero(cell(3,1)) .and. & 775 is_zero(cell(3,2)) ) then 776 cell_angle(3) = .true. 777 else 778 if ( IONOde ) then 779 write(*,'(a)') 'fixed: gamma-angle cannot & 780 &be constrained due to A and B having & 781 &components along Cartesian Z-direction.' 782 end if 783 call die('fixed: gamma-angle cannot be constrained.') 784 end if 785 786 end if 787 end do 788 789 else if ( leqi(namec,'cellside') .or. & 790 leqi(namec,'cell-side') .or. & 791 leqi(namec,'cellvector') .or. & 792 leqi(namec,'cell-vector') ) then 793 794 N = fdf_bnnames(pline) 795 796 do i = 2 , N 797 namec = fdf_bnames(pline,i) 798 ! Fix alpha angle 799 if ( leqi(namec,'a1') .or. leqi(namec,'a') ) then 800 if ( is_zero(cell(2,1)) .and. is_zero(cell(3,1)) ) then 801 cell_vector(1) = .true. 802 else 803 if ( IONOde ) then 804 write(*,'(a)') 'fixed: A-vector cannot & 805 &be constrained due to A having Y and Z & 806 &components.' 807 end if 808 call die('fixed: A-vector cannot be constrained.') 809 end if 810 811 else if ( leqi(namec,'a2') .or. leqi(namec,'b') ) then 812 if ( is_zero(cell(1,2)) .and. is_zero(cell(3,2)) ) then 813 cell_vector(2) = .true. 814 else 815 if ( IONOde ) then 816 write(*,'(a)') 'fixed: B-vector cannot & 817 &be constrained due to B having X and Z & 818 &components.' 819 end if 820 call die('fixed: B-vector cannot be constrained.') 821 end if 822 823 else if ( leqi(namec,'a3') .or. leqi(namec,'c') ) then 824 if ( is_zero(cell(1,3)) .and. is_zero(cell(2,3)) ) then 825 cell_vector(3) = .true. 826 else 827 if ( IONOde ) then 828 write(*,'(a)') 'fixed: C-vector cannot & 829 &be constrained due to C having X and Y & 830 &components.' 831 end if 832 call die('fixed: C-vector cannot be constrained.') 833 end if 834 835 end if 836 end do 837 838 ! ****** Now we only look at atomic specifications for constraints ****** 839 840 else if ( leqi(namec,'clear') .or. & 841 leqi(namec,'clear-previous') .or. leqi(namec,'clear-prev') ) then 842 843 ! This is a special option which removes all atoms in 844 ! the list for all currently found constraints (or the previous). 845 ! This lets one do pretty complex constraints. 846 847 ! For instance if one have a system of Au -- C -- Au with a 848 ! C being a hexagon attached to one Au on both sides, then bulk Au. 849 ! Then one can relax that one Au and the C by this constraint 850 ! Z 79 851 ! clear 11 45 852 call fix_brange(pline,rr,na) 853 854 ! Loop through all fixes that has been requested 855 sfix = 1 856 if ( leqi(namec,'clear-previous') ) sfix = ifix 857 if ( leqi(namec,'clear-prev') ) sfix = ifix 858 do ix = sfix , ifix 859 ! easy skip (easier than if) 860 if ( rr%n == 0 ) cycle 861 862 ! Remove all cleared atoms 863 call rgn_list(r_tmp,fixs(ix)%n,fixs(ix)%a) 864 call rgn_complement(rr,r_tmp,r_tmp) 865 if ( r_tmp%n /= fixs(ix)%n ) then 866 deallocate(fixs(ix)%a) 867 fixs(ix)%n = r_tmp%n 868 allocate(fixs(ix)%a(r_tmp%n)) 869 fixs(ix)%a(:) = r_tmp%r(:) 870 end if 871 call rgn_delete(r_tmp) 872 873 end do 874 875 else if ( leqi(namec,'position') .or. leqi(namec,'atom') ) then 876 877 ifix = ifix + 1 878 879 ! Create a list of atoms from this line 880 call fix_brange(pline,rr,na) 881 882 fixs(ifix)%n = rr%n 883 allocate(fixs(ifix)%a(rr%n)) 884 fixs(ifix)%a(:) = rr%r(:) 885 886 add_dir = .true. 887 fixs(ifix)%type = 'pos' 888 889 else if ( leqi(namec,'species-i') .or. leqi(namec,'Z') ) then 890 891 ifix = ifix + 1 892 893 ! Create a list of atoms from this line 894 ! We truncate to 1000 895 ! I.e. maximally use 1000 different species 896 ! Note that fix_brange also removes duplicates 897 call fix_brange(pline,rr,1000) 898 899 ! Loop through and allocate the correct atoms 900 ix = 0 901 N = 0 902 if ( leqi(namec,'species-i') ) then 903 do i = 1 , rr%n 904 N = N + count( isa(:) == rr%r(i) ) 905 end do 906 fixs(ifix)%n = N 907 allocate(fixs(ifix)%a(N)) 908 909 do i = 1 , na 910 if ( any( isa(i) == rr%r ) ) then 911 ix = ix + 1 912 fixs(ifix)%a(ix) = i 913 end if 914 end do 915 916 else if ( leqi(namec,'Z') ) then 917 do i = 1 , rr%n 918 N = N + count( iza(:) == rr%r(i) ) 919 end do 920 fixs(ifix)%n = N 921 allocate(fixs(ifix)%a(N)) 922 923 do i = 1 , na 924 if ( any( iza(i) == rr%r ) ) then 925 ix = ix + 1 926 fixs(ifix)%a(ix) = i 927 end if 928 end do 929 930 end if 931 932 add_dir = .true. 933 fixs(ifix)%type = 'pos' 934 935 else if ( leqi(namec,'rigid') .or. leqi(namec,'molecule') ) then 936 937 ! We restrict the entire molecule to move "together" 938 939 ifix = ifix + 1 940 941 ! Create a list of atoms from this line 942 call fix_brange(pline,rr,na) 943 944 fixs(ifix)%n = rr%n 945 allocate(fixs(ifix)%a(rr%n)) 946 fixs(ifix)%a(:) = rr%r(:) 947 948 add_dir = .true. 949 fixs(ifix)%type = 'rigid' 950 951 else if ( leqi(namec,'rigid-max') .or. leqi(namec,'molecule-max') ) then 952 953 ! We restrict the entire molecule to move "together" 954 955 ifix = ifix + 1 956 957 ! Create a list of atoms from this line 958 call fix_brange(pline,rr,na) 959 960 fixs(ifix)%n = rr%n 961 allocate(fixs(ifix)%a(rr%n)) 962 fixs(ifix)%a(:) = rr%r(:) 963 964 add_dir = .true. 965 fixs(ifix)%type = 'rigid-max' 966 967 else if ( leqi(namec,'center-of-mass') ) then 968 969 ! We restrict the center-of-mass for the specified region 970 ! of atoms 971 972 ifix = ifix + 1 973 974 ! Create a list of atoms from this line 975 call fix_brange(pline,rr,na) 976 977 fixs(ifix)%n = rr%n 978 allocate(fixs(ifix)%a(rr%n)) 979 fixs(ifix)%a(:) = rr%r(:) 980 981 fixs(ifix)%type = 'com' 982 983 else if ( leqi(namec,'center') ) then 984 985 ! We restrict the entire molecule to have the same center 986 ! of coordinates 987 988 ifix = ifix + 1 989 990 ! Create a list of atoms from this line 991 if ( (fdf_bntokens(pline) == 4 .and. fdf_bnreals(pline) == 3) & 992 .or. fdf_bntokens(pline) == 1 ) then 993 ! this line specifies all atoms (shorthand) 994 call rgn_range(rr,1,na) 995 else 996 call fix_brange(pline,rr,na) 997 end if 998 999 fixs(ifix)%n = rr%n 1000 allocate(fixs(ifix)%a(rr%n)) 1001 fixs(ifix)%a(:) = rr%r(:) 1002 1003 add_dir = .true. 1004 fixs(ifix)%type = 'center' 1005 1006 else if ( IONode ) then 1007 1008 write(*,'(2a)') 'WARNING: Skipping unrecognized geometry & 1009 &constraint: ',trim(namec) 1010 1011 end if 1012 1013 if ( add_dir ) then 1014 N = fdf_bnreals(pline) 1015 if ( N == 3 ) then 1016 1017 ! Fix the direction specified by the 3 reals 1018 fixs(ifix)%type = trim(fixs(ifix)%type)//'-dir' 1019 1020 fixs(ifix)%fix(1) = fdf_breals(pline,1) 1021 fixs(ifix)%fix(2) = fdf_breals(pline,2) 1022 fixs(ifix)%fix(3) = fdf_breals(pline,3) 1023 1024 fixs(ifix)%fix(:) = fixs(ifix)%fix(:) / VNORM( fixs(ifix)%fix(:) ) 1025 1026 else if ( N /= 0 ) then 1027 1028 write(*,'(a,i0,a)')'ERROR: Constraint reading: ',ifix,' was erroneous.' 1029 call die('You *must* specify 0 or 3 real values (not integers) & 1030 &to do a constraint in certain directions.') 1031 1032 end if 1033 end if 1034 1035 call rgn_delete(rr) 1036 1037 end do 1038 1039#ifdef DEBUG 1040 call write_debug( ' POS init_fixed' ) 1041#endif 1042 1043 contains 1044 1045 function is_ortho(v1,v2) result(is) 1046 real(dp), intent(in) :: v1(3), v2(3) 1047 logical :: is 1048 1049 is = is_zero( vnorm( vec_proj(v1,v2) ) ) 1050 1051 end function is_ortho 1052 1053 function is_zero(v) result(is) 1054 real(dp), intent(in) :: v 1055 logical :: is 1056 1057 is = abs(v) < 1.e-9_dp 1058 1059 end function is_zero 1060 1061 subroutine fix_brange(pline,r,na) 1062 type(parsed_line), pointer :: pline 1063 type(tRgn), intent(inout) :: r 1064 integer, intent(in) :: na 1065 integer :: i 1066 character(len=20) :: opt 1067 1068 ! If the number of names in this region is more than 1 1069 ! it could be an "all" designation 1070 if ( fdf_bnnames(pline) >= 2 ) then 1071 opt = fdf_bnames(pline,2) 1072 else 1073 opt = 'none' 1074 end if 1075 1076 if ( leqi(opt,'all') ) then 1077 1078 ! The user range is everything... 1079 call rgn_range(r,1,na) 1080 return 1081 1082 end if 1083 1084 ! Get entire range ( from -na to 2 * na ) 1085 ! we allow both negative and positive wrap-arounds 1086 ! (limit of one wrap-around) 1087 call fdf_brange(pline,r, 1, na) 1088 1089 ! This removes duplicates, then sorts it 1090 call rgn_uniq(r) 1091 call rgn_sort(r) 1092 1093 ! Remove zero 1094 i = rgn_pivot(r, 0) 1095 if ( i > 0 ) then 1096 i = rgn_pop(r, idx = i) 1097 end if 1098 1099 end subroutine fix_brange 1100 1101 end subroutine init_fixed 1102 1103 function is_fixed(ia) result(fixed) 1104 integer, intent(in) :: ia 1105 logical :: fixed, fi(3) 1106 integer :: if, i, N 1107 1108 fixed = .false. 1109 fi(:) = .false. 1110 1111 do if = 1 , N_fix 1112 1113 N = fixs(if)%n 1114 1115 do i = 1 , N 1116 1117 if ( ia == fixs(if)%a(i) ) then 1118 1119 if ( fixs(if)%type == 'pos' ) then 1120 fi(:) = .true. 1121 else if ( fixs(if)%fix(1) == 1._dp ) then 1122 fi(1) = .true. 1123 else if ( fixs(if)%fix(2) == 1._dp ) then 1124 fi(2) = .true. 1125 else if ( fixs(if)%fix(3) == 1._dp ) then 1126 fi(3) = .true. 1127 end if 1128 1129 exit 1130 1131 end if 1132 1133 end do 1134 1135 fixed = all(fi) 1136 if ( fixed ) return 1137 1138 end do 1139 1140 end function is_fixed 1141 1142 function is_constr(ia,type) result(constr) 1143 integer, intent(in) :: ia 1144 character(len=*), intent(in), optional :: type 1145 logical :: constr 1146 integer :: if, i, N 1147 1148 constr = .false. 1149 1150 do if = 1 , N_fix 1151 1152 N = fixs(if)%n 1153 1154 do i = 1 , N 1155 1156 if ( ia == fixs(if)%a(i) ) then 1157 if ( present(type) ) then 1158 constr = fixs(if)%type == type 1159 else 1160 constr = .true. 1161 end if 1162 if ( constr ) return 1163 end if 1164 1165 end do 1166 1167 end do 1168 1169 end function is_constr 1170 1171end module m_fixed 1172