1! This file is part of xtb. 2! 3! Copyright (C) 2017-2020 Stefan Grimme 4! 5! xtb is free software: you can redistribute it and/or modify it under 6! the terms of the GNU Lesser General Public License as published by 7! the Free Software Foundation, either version 3 of the License, or 8! (at your option) any later version. 9! 10! xtb is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU Lesser General Public License for more details. 14! 15! You should have received a copy of the GNU Lesser General Public License 16! along with xtb. If not, see <https://www.gnu.org/licenses/>. 17 18!! ======================================================================== 19! * WELCOME TO THE C O N S T R A I N S & S C A N S MODULE IN XTB * 20!! ------------------------------------------------------------------------ 21! The syntax is: 22!> $constrain 23!> ## this part is read by setparam.f90 24!> force constant=<real> 25!> all bonds=<bool> 26!> all angles=<bool> 27!> all torsions=<bool> 28!> ## this part is read in this module 29!> distance: <i>,<j>, auto 30!> distance: <i>,<j>, <real> 31!> angle: <i>,<j>,<k>, auto 32!> angle: <i>,<j>,<k>, <real> 33!> dihedral: <i>,<j>,<k>,<l>,auto 34!> dihedral: <i>,<j>,<k>,<l>,<real> 35!> center: <int>,<real> 36!> $fix 37!> force constant=<real> 38!> spring exponent=<int> 39!> atoms: <list> 40!> $split 41!> fragment1: <list> 42!> fragment2: <list> 43!> $wall 44!> sphere: auto,all 45!> sphere: auto,<list> 46!> sphere: <real>,all 47!> sphere: <real>,<list> 48!> ellipsoid: auto,all 49!> ellipsoid: auto,<list> 50!> ellipsoid: <real>,<real>,<real>,all 51!> ellipsoid: <real>,<real>,<real>,<list> 52!> $scan 53!> ... 54!> $end 55!! ======================================================================== 56module xtb_constrain_param 57 use xtb_mctc_accuracy, only : wp 58 use xtb_mctc_strings, only : parse 59 use xtb_readin, only : getline => strip_line,getValue,getListValue 60 use xtb_setparam, only : verbose 61 use xtb_type_environment, only : TEnvironment 62 use xtb_type_identitymap, only : TIdentityMap, init 63 use xtb_type_molecule, only : TMolecule 64 65 implicit none 66 67 private :: getline,getValue, handlerInterface 68 69 character,private,parameter :: flag = '$' 70 character,private,parameter :: colon = ':' 71 character,private,parameter :: space = ' ' 72 character,private,parameter :: equal = '=' 73 character,private,parameter :: hash = '#' 74 character,private,parameter :: comma = ',' 75 character,private,parameter :: semicolon = ';' 76 character,private,parameter :: dot = '.' 77 character(len=*),private,parameter :: flag_end = flag//'end' 78 79! Using allocatable arrays of dynamic length strings is only possible 80! with a lot of hacks, so we use good'ol fixed size stack arrays. 81! Let's choose something different from 42 that is not dividable by 10... ;) 82! Happy debugging! 83 integer,private,parameter :: p_str_length = 48 84 integer,private,parameter :: p_arg_length = 24 85 86 public 87 88 abstract interface 89 subroutine handlerInterface(env,key,val,nat,at,idMap,xyz) 90 import :: wp, TEnvironment, TIdentityMap 91 type(TEnvironment), intent(inout) :: env 92 character(len=*),intent(in) :: key 93 character(len=*),intent(in) :: val 94 integer, intent(in) :: nat 95 type(TIdentityMap), intent(in) :: idMap 96 integer, intent(in) :: at(nat) 97 real(wp),intent(in) :: xyz(3,nat) 98 end subroutine handlerInterface 99 end interface 100 101contains 102 103subroutine read_userdata(fname,env,mol) 104 use xtb_readin, only : find_new_name 105 use xtb_scanparam 106 implicit none 107 character(len=*), parameter :: source = 'userdata_read' 108 type(TEnvironment), intent(inout) :: env 109 type(TMolecule), intent(inout) :: mol 110 character(len=*),intent(in) :: fname 111 character(len=:),allocatable :: line 112 character(len=:),allocatable :: key 113 character(len=:),allocatable :: val 114 character(len=:),allocatable :: newname 115 type(TIdentityMap) :: idMap 116 integer :: i 117 integer :: id 118 integer :: ic 119 integer :: ie 120 integer :: err 121 logical :: exist 122 123 if (verbose) then 124 write(env%unit,'(72("$"))') 125 write(env%unit,'(1x,"CONSTRAINTS & SCANS: DEBUG SECTION")') 126 write(env%unit,'(72("$"))') 127 endif 128 129 call open_file(id,fname,'r') 130 if (id.eq.-1) then 131 call env%warning("could not find '"//fname//"'",source) 132 return 133 endif 134 rewind(id) ! not sure if this is necessary 135 136 call init(idMap, mol) 137 138 call idMap%writeInfo(env%unit) 139 140! read first line before the readloop starts, I have to do this 141! to avoid using backspace on id (dammit Turbomole format) 142 call getline(id,line,err) 143 readflags: do 144 ! check if there is a $ in the *first* column 145 if (index(line,flag).eq.1) then 146 select case(line(2:)) 147 case('fix' ) 148 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 149 call rdblock(env,set_fix, line,id,mol%n,mol%at,idMap,mol%xyz,err) 150 case('split' ) 151 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 152 call rdblock(env,set_split, line,id,mol%n,mol%at,idMap,mol%xyz,err) 153 case('constrain') 154 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 155 if (allocated(potset%xyz)) then 156 call rdblock(env,set_constr, line,id,mol%n,mol%at,idMap,potset%xyz,err) 157 else 158 call rdblock(env,set_constr, line,id,mol%n,mol%at,idMap,mol%xyz,err) 159 endif 160 case('scan' ) 161 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 162 call rdblock(env,set_scan, line,id,mol%n,mol%at,idMap,mol%xyz,err) 163 case('wall' ) 164 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 165 call rdblock(env,set_wall, line,id,mol%n,mol%at,idMap,mol%xyz,err) 166 case('metadyn' ) 167 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 168 call rdblock(env,set_metadyn,line,id,mol%n,mol%at,idMap,mol%xyz,err) 169 case('hess' ) 170 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 171 call rdblock(env,set_hess, line,id,mol%n,mol%at,idMap,mol%xyz,err) 172 case('path' ) 173 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 174 call rdblock(env,set_path, line,id,mol%n,mol%at,idMap,mol%xyz,err) 175 case('reactor' ) 176 if (verbose) write(env%unit,'(">",1x,a)') line(2:) 177 call rdblock(env,set_reactor,line,id,mol%n,mol%at,idMap,mol%xyz,err) 178 case('set' ); call rdsetbl(env,set_legacy,line,id,mol%n,mol%at,idMap,mol%xyz,err) 179 case default ! unknown keyword -> ignore, we don't raise them 180 call getline(id,line,err) 181 end select 182 else ! not a keyword -> ignore 183 call getline(id,line,err) 184 endif 185 ! check for end of file, which I will tolerate as alternative to $end 186 if (is_iostat_end(err)) exit readflags 187! if (index(line,flag_end).ne.0) exit readflags ! compatibility reasons 188 enddo readflags 189 190 if (verbose) write(env%unit,'(72("$"))') 191 call close_file(id) 192end subroutine read_userdata 193 194subroutine rdsetbl(env,handler,line,id,nat,at,idMap,xyz,err) 195 implicit none 196 character(len=*), parameter :: source = 'userdata_rdsetbl' 197 type(TEnvironment), intent(inout) :: env 198 integer,intent(in) :: id 199 procedure(handlerInterface) :: handler 200 integer, intent(in) :: nat 201 integer, intent(in) :: at(nat) 202 type(TIdentityMap), intent(in) :: idMap 203 real(wp),intent(in) :: xyz(3,nat) 204 integer,intent(out) :: err 205 character(len=:),allocatable,intent(out) :: line 206 character(len=:),allocatable :: key 207 character(len=:),allocatable :: val 208 integer :: ie 209 logical :: exitRun 210 do 211 call getline(id,line,err) 212 if (is_iostat_end(err)) exit 213 if (index(line,flag).ne.0) exit 214 if (verbose) write(env%unit,'("->",1x,a)') line 215 216 ! find the first colon 217 ie = index(line,space) 218 if ((line.eq.'').or.(ie.eq.0)) cycle 219 key = trim(line(:ie-1)) 220 val = trim(adjustl(line(ie+1:))) 221 call handler(env,key,val,nat,at,idMap,xyz) 222 call env%check(exitRun) 223 if (exitRun) then 224 call env%error("handler could not process input", source) 225 return 226 end if 227 enddo 228 229end subroutine rdsetbl 230 231subroutine rdblock(env,handler,line,id,nat,at,idMap,xyz,err) 232 implicit none 233 character(len=*), parameter :: source = 'userdata_rdblock' 234 type(TEnvironment), intent(inout) :: env 235 integer,intent(in) :: id 236 procedure(handlerInterface) :: handler 237 integer, intent(in) :: nat 238 integer, intent(in) :: at(nat) 239 type(TIdentityMap), intent(in) :: idMap 240 real(wp),intent(in) :: xyz(3,nat) 241 integer,intent(out) :: err 242 character(len=:),allocatable,intent(out) :: line 243 character(len=:),allocatable :: key 244 character(len=:),allocatable :: val 245 integer :: ie 246 logical :: exitRun 247 do 248 call getline(id,line,err) 249 if (is_iostat_end(err)) exit 250 if (index(line,flag).ne.0) exit 251 if (verbose) write(env%unit,'("->",1x,a)') line 252 253 ! find the first colon 254 ie = index(line,colon) 255 if ((line.eq.'').or.(ie.eq.0)) cycle 256 key = trim(line(:ie-1)) 257 val = trim(adjustl(line(ie+1:))) 258 call handler(env,key,val,nat,at,idMap,xyz) 259 call env%check(exitRun) 260 if (exitRun) then 261 call env%error("handler could not process input", source) 262 return 263 end if 264 enddo 265 266end subroutine rdblock 267 268subroutine set_fix(env,key,val,nat,at,idMap,xyz) 269 use xtb_type_atomlist 270 use xtb_fixparam 271 use xtb_setparam 272 implicit none 273 character(len=*), parameter :: source = 'userdata_fix' 274 type(TEnvironment), intent(inout) :: env 275 character(len=*),intent(in) :: key 276 character(len=*),intent(in) :: val 277 integer, intent(in) :: nat 278 integer, intent(in) :: at(nat) 279 type(TIdentityMap), intent(in) :: idMap 280 real(wp),intent(in) :: xyz(3,nat) 281 282 type(TAtomList) :: atl 283 284 integer :: i 285 integer :: iat 286 integer :: idum 287 integer :: nlist 288 integer, allocatable :: list(:) 289 real(wp) :: ddum 290 logical :: ldum 291 292 integer :: narg 293 character(len=p_str_length),dimension(p_arg_length) :: argv 294 295 call atl%resize(nat) 296 297 call parse(val,comma,argv,narg) 298! some debug xtb_printout 299 if (verbose) then 300 do idum = 1, narg 301 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 302 enddo 303 endif 304 select case(key) 305 case default ! ignore, don't even think about raising them 306 case('elements') 307 call atl%new 308 do idum = 1, narg 309 ! get element by symbol 310 if (idMap%has(argv(idum))) then 311 call idMap%get(list, argv(idum)) 312 if (allocated(list)) then 313 call atl%add(list) 314 else 315 call env%warning("Unknown element: '"//trim(argv(idum))//"'",source) 316 cycle 317 end if 318 else 319 ldum = getValue(env,trim(argv(idum)),iat) 320 if (.not.ldum) cycle ! skip garbage input 321 ! check for unreasonable input 322 if (iat > 0) then 323 call atl%add(at.eq.iat) 324 else 325 call env%warning("Unknown element: '"//trim(argv(idum))//"'",source) 326 cycle 327 endif 328 endif 329 enddo 330 if (fixset%n > 0) call atl%add(fixset%atoms(:fixset%n)) 331 call atl%to_list(list) 332 fixset%atoms = list 333 fixset%n = size(list) 334 case('atoms') 335 call atl%new(val) 336 if (atl%get_error()) then 337 call env%warning('something is wrong in the fixing list',source) 338 return 339 endif 340 if (fixset%n > 0) call atl%add(fixset%atoms(:fixset%n)) 341 call atl%to_list(list) 342 fixset%atoms = list 343 fixset%n = size(list) 344 case('freeze') 345 call atl%new(val) 346 if (atl%get_error()) then 347 call env%warning('something is wrong in the freezing list',source) 348 return 349 endif 350 if (freezeset%n > 0) call atl%add(freezeset%atoms(:freezeset%n)) 351 call atl%to_list(list) 352 freezeset%atoms = list 353 freezeset%n = size(list) 354 case('shake') 355 allocate(list(nat*(nat+1)/2), source=0) 356 if (mod(narg,2).ne.0) then 357 call env%warning("could not read input for user defined shake!",source) 358 return 359 endif 360 if (narg+shakeset%n > nat*(nat+1)/2) then 361 call env%warning("too many SHAKE constraints!",source) 362 return 363 endif 364 if (.not.shake_md) shake_md = .true. 365 do idum = 1, narg 366 if (getValue(env,trim(argv(idum)),iat)) then 367 if (iat.gt.nat) then 368 call env%warning('Attempted constrain atom not present in molecule.',source) 369 cycle 370 endif 371 shakeset%n = shakeset%n+1 372 shakeset%atoms(shakeset%n) = iat 373 else 374 call env%warning("Something went wrong in set_fix_ 'shake'.",source) 375 return ! you screwed it, let's get out of here 376 endif 377 enddo 378 end select 379 380end subroutine set_fix 381 382subroutine set_constr(env,key,val,nat,at,idMap,xyz) 383 use xtb_mctc_constants 384 use xtb_mctc_convert 385 use xtb_type_atomlist 386 use xtb_scanparam 387 use xtb_splitparam 388 implicit none 389 character(len=*), parameter :: source = 'userdata_constr' 390 type(TEnvironment), intent(inout) :: env 391 character(len=*),intent(in) :: key 392 character(len=*),intent(in) :: val 393 integer, intent(in) :: nat 394 integer, intent(in) :: at(nat) 395 type(TIdentityMap), intent(in) :: idMap 396 real(wp),intent(in) :: xyz(3,nat) 397 398 type(TAtomList) :: atl 399 400 integer :: iat 401 integer :: ioffset 402 integer :: idum 403 real(wp) :: ddum 404 integer :: nlist 405 integer, allocatable :: list(:) 406 logical :: ldum 407 integer :: i,j,k,l 408 real(wp) :: phi,dist,ra(3),rb(3) 409 real(wp),external :: valijkl 410 411 integer :: narg 412 character(len=p_str_length),dimension(p_arg_length) :: argv 413 414 call atl%resize(nat) 415 416 call parse(val,comma,argv,narg) 417 if (verbose) then 418 do idum = 1, narg 419 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 420 enddo 421 endif 422 select case(key) 423 case default ! ignore, don't even think about raising them 424 425 case('elements') 426 call atl%new 427 do idum = 1, narg 428 ! get element by symbol 429 if (idMap%has(argv(idum))) then 430 call idMap%get(list, argv(idum)) 431 if (allocated(list)) then 432 call atl%add(list) 433 else 434 call env%warning("Unknown element: '"//trim(argv(idum))//"'",source) 435 cycle 436 end if 437 else 438 ldum = getValue(env,trim(argv(idum)),iat) 439 if (.not.ldum) cycle ! skip garbage input 440 ! check for unreasonable input 441 if (iat > 0) then 442 call atl%add(at.eq.iat) 443 else 444 call env%warning("Unknown element: '"//trim(argv(idum))//"'",source) 445 cycle 446 endif 447 endif 448 enddo 449 if (potset%pos%n > 0) call atl%add(potset%pos%atoms(:potset%pos%n)) 450 call atl%to_list(list) 451 potset%pos%atoms = list 452 potset%pos%n = size(list) 453 case('atoms') 454 call atl%new(val) 455 if (atl%get_error()) then 456 call env%warning('something is wrong in the fixing list',source) 457 return 458 endif 459 if (potset%pos%n > 0) call atl%add(potset%pos%atoms(:potset%pos%n)) 460 call atl%to_list(list) 461 potset%pos%atoms = list 462 potset%pos%n = size(list) 463 464 case('DISTANCE') 465 if (narg.lt.3 .or. narg.gt.4) then 466 call env%error('not enough arguments to constrain a distance',source) 467 return 468 endif 469 ioffset = 2*potset%dist%n 470 potset%dist%n = potset%dist%n+1 471 ! part 1: get the constrained atoms 472 do i = 1, 2 473 if (getValue(env,trim(argv(i)),idum)) then 474 potset%dist%atoms(ioffset+i) = idum 475 else 476 call env%warning("Something went wrong in set_constr_ 'distance'. (1)",source) 477 potset%dist%n = potset%dist%n-1 ! remove invalid contrain 478 return 479 endif 480 enddo 481 ! part 2: get the distance between those atoms 482 i = potset%dist%atoms(ioffset+1) 483 j = potset%dist%atoms(ioffset+2) 484 dist = sqrt(sum((xyz(:,i)-xyz(:,j))**2)) 485 if (trim(argv(3)).eq.'auto') then 486 potset%dist%val(potset%dist%n) = dist 487 else 488 if (getValue(env,trim(argv(3)),ddum)) then 489 potset%dist%val(potset%dist%n) = ddum * aatoau 490 else 491 call env%warning("Something went wrong in set_constr_ 'distance'. (2)",source) 492 potset%dist%n = potset%dist%n-1 ! remove invalid contrain 493 return 494 endif 495 endif 496 if (narg.eq.4) then 497 if (getValue(env,trim(argv(4)),idum)) then 498 if (idum < 2 .or. mod(idum, 2).ne.0) then 499 call env%warning("Invalid spring exponent given", source) 500 potset%dist%n = potset%dist%n-1 ! remove invalid contrain 501 return 502 end if 503 potset%dist%expo(potset%dist%n) = real(idum, wp) 504 else 505 call env%warning("Something went wrong in set_constr_ 'distance'. (3)",source) 506 potset%dist%n = potset%dist%n-1 ! remove invalid contrain 507 return 508 end if 509 write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//& 510 '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å",1x,"with expo",1x,i0)') & 511 i,j, potset%dist%val(potset%dist%n)*autoaa, dist*autoaa, & 512 nint(potset%dist%expo(potset%dist%n)) 513 else 514 write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//& 515 '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å")') & 516 i,j, potset%dist%val(potset%dist%n)*autoaa, dist*autoaa 517 end if 518 519 case('ANGLE') 520 if (narg.ne.4) then 521 call env%error('not enough arguments to constrain an angle',source) 522 return 523 endif 524 ioffset = 3*potset%angle%n 525 potset%angle%n = potset%angle%n+1 526 ! part 1: get the constrained atoms 527 do i = 1, narg-1 528 if (getValue(env,trim(argv(i)),idum)) then 529 potset%angle%atoms(ioffset+i) = idum 530 else 531 call env%warning("Something went wrong in set_constr_ 'angle'. (1)",source) 532 potset%angle%n = potset%angle%n-1 ! remove invalid contrain 533 return 534 endif 535 enddo 536 ! part 2: get the angle between the constrained atoms 537 i = potset%angle%atoms(ioffset+1) 538 j = potset%angle%atoms(ioffset+2) 539 k = potset%angle%atoms(ioffset+3) 540 call bangl(xyz,i,j,k,phi) 541 if (trim(argv(narg)).eq.'auto') then 542 potset%angle%val(potset%angle%n) = phi 543 else 544 if (getValue(env,trim(argv(narg)),ddum)) then 545 potset%angle%val(potset%angle%n) = pi/180.0_wp * ddum 546 else 547 call env%warning("Something went wrong in set_constr_ 'angle'. (2)",source) 548 potset%angle%n = potset%angle%n-1 ! remove invalid contrain 549 return 550 endif 551 endif 552 553 case('DIHEDRAL') 554 if (narg.ne.5) then 555 call env%error('not enough arguments to constrain a dihedral',source) 556 return 557 endif 558 ioffset = 4*potset%dihedral%n 559 potset%dihedral%n = potset%dihedral%n+1 560 if (nconstr.gt.maxconstr) then ! double check this 561 call env%error('This should never happen! Let somebody check set_constr',source) 562 return 563 end if 564 ! part 1: get the constrained atoms 565 do i = 1, narg-1 566 if (getValue(env,trim(argv(i)),idum)) then 567 potset%dihedral%atoms(ioffset+i) = idum 568 else 569 call env%warning("Something went wrong in set_constr_ 'dihedral'. (1)",source) 570 potset%dihedral%n = potset%dihedral%n-1 ! remove invalid contrain 571 return 572 endif 573 enddo 574 ! part 2: get the angle between the constrained atoms 575 i = potset%dihedral%atoms(ioffset+1) 576 j = potset%dihedral%atoms(ioffset+2) 577 k = potset%dihedral%atoms(ioffset+3) 578 l = potset%dihedral%atoms(ioffset+4) 579 phi=valijkl(nat,xyz,i,j,k,l) 580 if (trim(argv(narg)).eq.'auto') then 581 potset%dihedral%val(potset%dihedral%n) = phi 582 else 583 if (getValue(env,trim(argv(narg)),ddum)) then 584 potset%dihedral%val(potset%dihedral%n) = pi/180.0_wp * ddum 585 else 586 call env%warning("Something went wrong in set_constr_ 'dihedral'. (2)",source) 587 potset%dihedral%n = potset%dihedral%n-1 ! remove invalid contrain 588 return 589 endif 590 endif 591 592 case('distance') 593 if (narg.ne.3) then 594 call env%error('not enough arguments to constrain a distance',source) 595 return 596 endif 597 nconstr = nconstr+1 598 if (nconstr.gt.maxconstr) then ! double check this 599 call env%error('This should never happen! Let somebody check set_constr',source) 600 return 601 end if 602 ! part 1: get the constrained atoms 603 do i = 1, narg-1 604 if (getValue(env,trim(argv(i)),idum)) then 605 atconstr(i,nconstr) = idum 606 else 607 call env%warning("Something went wrong in set_constr_ 'distance'. (1)",source) 608 nconstr = nconstr-1 ! remove invalid contrain 609 return 610 endif 611 enddo 612 ! part 2: get the distance between those atoms 613 i = atconstr(1,nconstr) 614 j = atconstr(2,nconstr) 615 dist = sqrt(sum((xyz(:,i)-xyz(:,j))**2)) 616 if (trim(argv(narg)).eq.'auto') then 617 valconstr(nconstr) = dist 618 else 619 if (getValue(env,trim(argv(narg)),ddum)) then 620 valconstr(nconstr) = ddum * aatoau 621 else 622 call env%warning("Something went wrong in set_constr_ 'distance'. (2)",source) 623 nconstr = nconstr-1 ! remove invalid contrain 624 return 625 endif 626 endif 627 write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//& 628 '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å")') & 629 i,j, valconstr(nconstr)*autoaa, dist*autoaa 630 631 case('angle') 632 if (narg.ne.4) then 633 call env%error('not enough arguments to constrain an angle',source) 634 return 635 endif 636 nconstr = nconstr+1 637 if (nconstr.gt.maxconstr) then ! double check this 638 call env%error('This should never happen! Let somebody check set_constr',source) 639 return 640 end if 641 ! part 1: get the constrained atoms 642 do i = 1, narg-1 643 if (getValue(env,trim(argv(i)),idum)) then 644 atconstr(i,nconstr) = idum 645 else 646 call env%warning("Something went wrong in set_constr_ 'angle'. (1)",source) 647 nconstr = nconstr-1 ! remove invalid contrain 648 return 649 endif 650 enddo 651 ! part 2: get the angle between the constrained atoms 652 i = atconstr(1,nconstr) 653 j = atconstr(2,nconstr) 654 k = atconstr(3,nconstr) 655 call bangl(xyz,i,j,k,phi) 656 if (trim(argv(narg)).eq.'auto') then 657 valconstr(nconstr) = phi 658 else 659 if (getValue(env,trim(argv(narg)),ddum)) then 660 valconstr(nconstr) = pi/180.0_wp * ddum 661 else 662 call env%warning("Something went wrong in set_constr_ 'angle'. (2)",source) 663 nconstr = nconstr-1 ! remove invalid contrain 664 return 665 endif 666 endif 667 write(env%unit,'("constraining angle",3(1x,i0),1x,"to",'//& 668 '1x,f12.7,"°, actual value:",1x,f12.7,"°")') & 669 i,j,k,180.0_wp/pi * valconstr(nconstr),180.0_wp/pi * phi 670 671 case('dihedral') 672 if (narg.ne.5) then 673 call env%error('not enough arguments to constrain a dihedral',source) 674 return 675 endif 676 nconstr = nconstr+1 677 if (nconstr.gt.maxconstr) then ! double check this 678 call env%error('This should never happen! Let somebody check set_constr',source) 679 return 680 end if 681 ! part 1: get the constrained atoms 682 do i = 1, narg-1 683 if (getValue(env,trim(argv(i)),idum)) then 684 atconstr(i,nconstr) = idum 685 else 686 call env%warning("Something went wrong in set_constr_ 'dihedral'. (1)",source) 687 nconstr = nconstr-1 ! remove invalid contrain 688 return 689 endif 690 enddo 691 ! part 2: get the angle between the constrained atoms 692 i = atconstr(1,nconstr) 693 j = atconstr(2,nconstr) 694 k = atconstr(3,nconstr) 695 l = atconstr(4,nconstr) 696 phi=valijkl(nat,xyz,i,j,k,l) 697 if (trim(argv(narg)).eq.'auto') then 698 valconstr(nconstr) = phi 699 else 700 if (getValue(env,trim(argv(narg)),ddum)) then 701 valconstr(nconstr) = pi/180.0_wp * ddum 702 else 703 call env%warning("Something went wrong in set_constr_ 'dihedral'. (2)",source) 704 nconstr = nconstr-1 ! remove invalid contrain 705 return 706 endif 707 endif 708 write(env%unit,'("constraining angle",4(1x,i0),1x,"to",'//& 709 '1x,f12.7,"°, actual value:",1x,f12.7,"°")') & 710 i,j,k,l,180.0_wp/pi * valconstr(nconstr),180.0_wp/pi * phi 711 712 case('center') 713 ! copied from molbld.f without modification (except for error handling) 714 if (narg.ne.2) then 715 call env%error('not enough argument to constrain center',source) 716 return 717 endif 718 nconstr = nconstr+1 719 atconstr(1,nconstr) = -2 720 if (getValue(env,trim(argv(2)),ddum)) then 721 valconstr(nconstr) = ddum 722 else 723 call env%warning("Something went wrong in set_constr_ 'center'. (1)",source) 724 nconstr = nconstr-1 725 return 726 endif 727 if (getValue(env,trim(argv(1)),idum)) then 728 iatf1 = idum 729 else 730 call env%warning("Something went wrong in set_constr_ 'center'. (2)",source) 731 nconstr = nconstr-1 732 return 733 endif 734 massf1 = 0.0_wp 735 do i = 1, iatf1 736 massf1 = massf1 + atmass(i) 737 enddo 738 739 case('cma','cma interface') 740 ! copied from molbld.f without modification (except for error handling) 741 if (narg.ne.1) then 742 call env%error('not enough argument to constrain cma',source) 743 return 744 endif 745 if (key.eq.'cma interface') call cmaiface(nat,at,xyz) 746 nconstr = nconstr+1 747 atconstr(1,nconstr) = 0 748 if (trim(argv(1)).eq.'auto') then 749 massf1 = 0.0_wp 750 massf2 = 0.0_wp 751 do i = 1, nat 752 if (splitlist(i).eq.1) then 753 massf1 = massf1 + atmass(i) 754 else 755 massf2 = massf2 + atmass(i) 756 endif 757 enddo 758 call cmafrag(nat,at,xyz,ra,rb) 759 valconstr(nconstr) = rcma 760 write(env%unit,'("constraining fragment CMA to initial R(Ang.)=",f8.3)')& 761 valconstr(nconstr)*autoaa 762 else 763 if (getValue(env,trim(argv(1)),ddum)) then 764 valconstr(nconstr) = ddum*aatoau 765 else 766 call env%warning("Something went wrong in set_constr_ 'cma'.",source) 767 nconstr = nconstr-1 768 return 769 endif 770 endif 771 772! case('x','y','z') 773! idum = index('xyz',key) 774 case('z') 775 ! copied from molbld.f without modification (except for error handling) 776 if (narg.ne.1) then 777 call env%error('not enough argument to constrain z coordinate',source) 778 return 779 endif 780 zconstr = 1 781 nconstr = nconstr+1 782 atconstr(1,nconstr) = -1 783 if (getValue(env,trim(argv(1)),ddum)) then 784 valconstr(nconstr) = ddum*aatoau 785 else 786 call env%warning("Something went wrong in set_constr_ 'z'.",source) 787 nconstr = nconstr-1 788 return 789 endif 790 if (sum(abs(xyz(3,1:iatf1))).gt.1.d-3) then 791 call env%warning('z-coordinates of fragment 1 must be 0!',source) 792 nconstr = nconstr-1 793 return 794 endif 795 write(env%unit,'("constraining fragment 2 to Z=0 plane at (Ang.)=",f8.3)') & 796 valconstr(nconstr)*autoaa 797 798 end select 799end subroutine set_constr 800 801!! --------------------------------------------------------------[SAW1809]- 802! this is the new version of the scan routine exploiting all features 803subroutine set_scan(env,key,val,nat,at,idMap,xyz) 804 use xtb_scanparam 805 use xtb_mctc_convert, only : aatoau 806 use xtb_mctc_constants, only : pi 807 implicit none 808 character(len=*), parameter :: source = 'userdata_scan' 809 type(TEnvironment), intent(inout) :: env 810 character(len=*),intent(in) :: key 811 character(len=*),intent(in) :: val 812 integer, intent(in) :: nat 813 integer, intent(in) :: at(nat) 814 type(TIdentityMap), intent(in) :: idMap 815 real(wp),intent(in) :: xyz(3,nat) 816 817 integer :: idum 818 real(wp) :: ddum 819 logical :: ldum 820 821 integer :: i,ie 822 real(wp) :: start_value,end_value 823 integer :: narg 824 character(len=p_str_length),dimension(p_arg_length) :: argv 825 character(len=:),allocatable :: temp 826 827 nscan = nscan+1 ! make a new scan 828 829 ie = index(val,semicolon) 830 if (ie.ne.0) then 831 temp = val(:ie-1) 832 idum = nconstr 833 call set_constr(env,key,temp,nat,at,idMap,xyz) ! generate a new constraint 834 scan_list(nscan)%iconstr = nconstr ! new generated 835 if (idum.eq.nconstr) then 836 call env%error('Failed to generate constraint',source) 837 return 838 end if 839 temp = val(1+ie:) 840 else 841 if (getValue(env,key,idum)) then 842 if (idum.gt.nconstr) then 843 call env%error('Constraint '''//key//''' is not defined',source) 844 return 845 endif 846 scan_list(nscan)%iconstr = idum 847 else 848 call env%error('Constraint '''//key//''' is invalid in this context',source) 849 return 850 endif 851 temp = val 852 endif 853 854 call parse(temp,comma,argv,narg) 855 if (verbose) then 856 do idum = 1, narg 857 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 858 enddo 859 endif 860 861 ! get back the constraint number which should be scanned 862 i = scan_list(nscan)%iconstr 863 ! instead of simply saving the kind of constraint while reading, 864 ! we used some really obvious integer code system depending on 865 ! the number of elements in the atconstr-Array. 866 ! Expected to break at some point, add a FIXME and wait for complaints... 867 idum = atconstr(1,i) 868 if (atconstr(2,i).gt.0) idum = 1 869 if (atconstr(3,i).gt.0) idum = 2 870 if (atconstr(4,i).gt.0) idum = 3 871 872 if (getValue(env,trim(argv(1)),ddum)) then 873 if (idum.le.1) then 874 start_value = ddum * aatoau 875 else 876 start_value = ddum * pi/180.0_wp 877 endif 878 else 879 call env%error('Invalid start value for scan',source) 880 return 881 endif 882 if (getValue(env,trim(argv(2)),ddum)) then 883 if (idum.le.1) then 884 end_value = ddum * aatoau 885 else 886 end_value = ddum * pi/180.0_wp 887 endif 888 else 889 call env%error('Invalid end value for scan',source) 890 return 891 endif 892 if (getValue(env,trim(argv(3)),idum)) then 893 scan_list(nscan)%nscan = idum 894 allocate( scan_list(nscan)%valscan(idum), source = 0.0_wp ) 895 else 896 call env%error('Invalid step number for scan',source) 897 return 898 endif 899 900 do i = 1, scan_list(nscan)%nscan 901 scan_list(nscan)%valscan(i) = start_value + (end_value-start_value) & 902 & * real(i-1,wp)/real(scan_list(nscan)%nscan-1,wp) 903 if (verbose) write(env%unit,'(i5,1x,f12.8)') i,scan_list(nscan)%valscan(i) 904 enddo 905 906end subroutine set_scan 907 908subroutine set_wall(env,key,val,nat,at,idMap,xyz) 909 use xtb_mctc_convert, only : autoaa 910 use xtb_sphereparam 911 implicit none 912 character(len=*), parameter :: source = 'userdata_wall' 913 type(TEnvironment), intent(inout) :: env 914 character(len=*),intent(in) :: key 915 character(len=*),intent(in) :: val 916 integer, intent(in) :: nat 917 integer, intent(in) :: at(nat) 918 type(TIdentityMap), intent(in) :: idMap 919 real(wp),intent(in) :: xyz(3,nat) 920 921 integer :: idum 922 real(wp) :: ddum,darray(3) 923 logical :: ldum 924 integer :: list(nat),nlist 925 integer :: tlist(nat),ntlist 926 integer :: iarg,iaxis 927 real(wp) :: radius,radii(3),center(3) 928 929 integer :: narg 930 character(len=p_str_length),dimension(p_arg_length) :: argv 931 932 nlist = 0 933 ntlist = 0 934 list = 0 935 tlist = 0 936 937 radius = 0.0_wp 938 radii = 0.0_wp 939 center = 0.0_wp 940 941 call parse(val,comma,argv,narg) 942 if (verbose) then 943 do idum = 1, narg 944 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 945 enddo 946 endif 947 select case(key) 948 case default ! ignore, don't even think about raising them 949 case('sphere') 950 if (narg.lt.2) then 951 call env%warning("Not enough arguments to set up a spherical wall",source) 952 return 953 endif 954 ! part 1: get the sphere radius 955 if (trim(argv(1)).eq.'auto') then 956 call get_sphere_radius(nat,at,xyz,center,radius,do_trafo=.false.) 957 else 958 if (getValue(env,trim(argv(1)),ddum)) then 959 radius = ddum 960 center = 0.0_wp 961 else 962 ! warning already generated by get_value 963 return ! something went wrong 964 endif 965 endif 966 ! part 2: get atoms 967 if (trim(argv(2)).eq.'all') then 968 call set_sphere_radius(radius,center) 969 else 970 do iarg = 2, narg 971 if (getListValue(env,trim(argv(iarg)),tlist,ntlist)) then 972 if (nlist+ntlist.gt.nat) then 973 call env%warning("Too many atoms in list for spherical wall.",source) 974 return ! something went wrong 975 endif 976 if (maxval(tlist(:ntlist)).gt.nat) then 977 call env%warning("Attempted to wall in a non-existing atom",source) 978 cycle ! skip crappy input 979 endif 980 list(nlist+1:nlist+ntlist) = tlist 981 nlist = nlist + ntlist 982 else 983 ! warning already generated by get_list_value 984 return ! something went wrong 985 endif 986 enddo 987 call set_sphere_radius(radius,center,nlist,list) 988 endif 989 write(env%unit,'("spherical wallpotenial with radius",'//& 990 '1x,f12.7,1x,"Å")') radius*autoaa 991 992 case('ellipsoid') 993 if (narg.lt.4) then 994 call env%warning("Not enough arguments to set up an ellipsoidal wall",source) 995 return 996 endif 997 ! part 1: get ellipsoid axis 998 if ((trim(argv(1)).eq.'auto').or.(trim(argv(2)).eq.'auto').or. & 999 (trim(argv(3)).eq.'auto')) then 1000 ! use get_ellipsoid_radius if you manage to implement it 1001 call get_sphere_radius(nat,at,xyz,center,radius,do_trafo=.false.) 1002 endif 1003 do iaxis = 1, 3 1004 if (trim(argv(iaxis)).eq.'auto') then 1005 radii(iaxis) = radius 1006 else 1007 if (getValue(env,trim(argv(iaxis)),ddum)) then 1008 radii(iaxis) = ddum 1009 center(iaxis) = 0.0_wp 1010 else 1011 ! warning already generated by get_value 1012 return ! something went wrong 1013 endif 1014 endif 1015 enddo 1016 ! part 2: get atoms 1017 if (trim(argv(4)).eq.'all') then 1018 call set_sphere_radius(radii,center) 1019 else 1020 do iarg = 4, narg 1021 if (getListValue(env,trim(argv(iarg)),tlist,ntlist)) then 1022 if (nlist+ntlist.gt.nat) then 1023 call env%warning("Too many atoms in list for spherical wall.",source) 1024 return ! something went wrong 1025 endif 1026 if (maxval(tlist(:ntlist)).gt.nat) then 1027 call env%warning("Attempted to wall in a non-existing atom",source) 1028 cycle ! skip crappy input 1029 endif 1030 list(nlist+1:nlist+ntlist) = tlist 1031 nlist = nlist + ntlist 1032 else 1033 ! warning already generated by get_list_value 1034 return ! something went wrong 1035 endif 1036 enddo 1037 call set_sphere_radius(radii,center,nlist,list) 1038 endif 1039 write(env%unit,'("ellipsoidal wallpotenial with radii",'//& 1040 '3(1x,f12.7,1x,"Å"))') radii*autoaa 1041 1042 end select 1043 1044end subroutine set_wall 1045 1046subroutine set_split(env,key,val,nat,at,idMap,xyz) 1047 use xtb_splitparam 1048 implicit none 1049 character(len=*), parameter :: source = 'userdata_split' 1050 type(TEnvironment), intent(inout) :: env 1051 character(len=*),intent(in) :: key 1052 character(len=*),intent(in) :: val 1053 integer, intent(in) :: nat 1054 integer, intent(in) :: at(nat) 1055 type(TIdentityMap), intent(in) :: idMap 1056 real(wp),intent(in) :: xyz(3,nat) 1057 1058 integer :: idum 1059 real(wp) :: ddum 1060 logical :: ldum 1061 integer :: list(nat),nlist 1062 integer :: i,j 1063 1064 integer :: narg 1065 character(len=p_str_length),dimension(p_arg_length) :: argv 1066 1067 call parse(val,comma,argv,narg) 1068 if (verbose) then 1069 do idum = 1, narg 1070 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 1071 enddo 1072 endif 1073 select case(key) 1074 case default ! ignore, don't even think about raising them 1075 case('fragment1') 1076 do i = 1, narg 1077 if (getListValue(env,trim(argv(i)),list,nlist)) then 1078 do j = 1, nlist 1079 splitlist(list(j)) = 1 1080 enddo 1081 iatf1 = iatf1 + nlist ! compatibility stuff, as usual 1082 iatf(1) = iatf(1) + nlist 1083 else 1084 return ! you screwed it, let's get out of here 1085 endif 1086 enddo 1087 case('fragment2') 1088 do i = 1, narg 1089 if (getListValue(env,trim(argv(i)),list,nlist)) then 1090 do j = 1, nlist 1091 splitlist(list(j)) = 2 1092 enddo 1093 iatf2 = iatf2 + nlist ! compatibility stuff, as usual 1094 iatf(2) = iatf(2) + nlist 1095 else 1096 return ! you screwed it, let's get out of here 1097 endif 1098 enddo 1099 case('fragment') 1100 if (getValue(env,trim(argv(1)),idum)) then 1101 if (idum.gt.nat) then 1102 call env%warning("rejecting fragment number greater than number of atoms",source) 1103 return ! doesn't really make sense, sorry, your problem 1104 endif 1105 do i = 2, narg 1106 if (getListValue(env,trim(argv(i)),list,nlist)) then 1107 do j = 1, nlist 1108 splitlist(list(j)) = idum 1109 enddo 1110 iatf(idum) = iatf(idum) + nlist 1111 if (idum.eq.1) then ! legacy, don't remove or something breaks 1112 iatf1 = iatf1 + nlist 1113 else if (idum.eq.2) then 1114 iatf2 = iatf2 + nlist 1115 endif 1116 else 1117 return ! you screwed it, let's get out of here 1118 endif 1119 enddo 1120 endif 1121 end select 1122end subroutine set_split 1123 1124subroutine set_hess(env,key,val,nat,at,idMap,xyz) 1125 use xtb_splitparam 1126 implicit none 1127 character(len=*), parameter :: source = 'userdata_hess' 1128 type(TEnvironment), intent(inout) :: env 1129 character(len=*),intent(in) :: key 1130 character(len=*),intent(in) :: val 1131 integer, intent(in) :: nat 1132 integer, intent(in) :: at(nat) 1133 type(TIdentityMap), intent(in) :: idMap 1134 real(wp),intent(in) :: xyz(3,nat) 1135 1136 integer :: idum 1137 real(wp) :: ddum 1138 logical :: ldum 1139 integer :: i,j 1140 integer, allocatable :: list(:) 1141 1142 integer :: narg 1143 character(len=p_str_length),dimension(p_arg_length) :: argv 1144 1145 call parse(val,comma,argv,narg) 1146 if (verbose) then 1147 do idum = 1, narg 1148 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 1149 enddo 1150 endif 1151 select case(key) 1152 case default ! ignore, don't even think about raising them 1153 case('element mass') 1154 if (mod(narg,2).ne.0) then 1155 call env%warning("Something went wrong in set_mass_ 'element'.",source) 1156 endif 1157 do i = 1, narg, 2 1158 j = i+1 1159 if (getValue(env,trim(argv(j)),ddum)) then 1160 if (idMap%has(argv(i))) then 1161 call idMap%set(atmass, argv(i), ddum) 1162 write(env%unit,'(a,a,a,1x,g0)') & 1163 "mass of elements '",trim(argv(i)),"' changed to", ddum 1164 else if (getValue(env,trim(argv(i)),idum)) then 1165 where(at.eq.idum) atmass = ddum 1166 write(env%unit,'(a,i0,a,1x,g0)') & 1167 "mass of elements with Z=",idum," changed to", ddum 1168 end if 1169 end if 1170 enddo 1171 case('modify mass','isotope') 1172 if (mod(narg,2).ne.0) then 1173 call env%warning("Something went wrong in set_mass_ 'modify'.",source) 1174 endif 1175 do i = 1, narg, 2 1176 j = i+1 1177 if (getValue(env,trim(argv(i)),idum).and.& 1178 getValue(env,trim(argv(j)),ddum)) then 1179 if (idum.gt.nat) then 1180 call env%warning('Attempted setting atom mass not present in system.',source) 1181 cycle 1182 endif 1183 atmass(idum) = ddum 1184 write(env%unit,'(a,1x,i0,1x,a,1x,g0)') & 1185 'mass of atom ',idum,' changed to',atmass(idum) 1186 endif 1187 enddo 1188 case('scale mass') 1189 if (mod(narg,2).ne.0) then 1190 call env%warning("Something went wrong in set_mass_ 'scale'.",source) 1191 endif 1192 do i = 1, narg, 2 1193 j = i+1 1194 if (getValue(env,trim(argv(i)),idum).and.& 1195 getValue(env,trim(argv(j)),ddum)) then 1196 if (idum.gt.nat) then 1197 call env%warning('Attempted scaling atom not present in system.',source) 1198 cycle 1199 endif 1200 atmass(idum) = atmass(idum)*ddum 1201 write(env%unit,'(a,1x,i0,1x,a,1x,g0)') & 1202 'mass of atom ',idum,' changed to',atmass(idum) 1203 endif 1204 enddo 1205 end select 1206 1207end subroutine set_hess 1208 1209subroutine set_reactor(env,key,val,nat,at,idMap,xyz) 1210 use xtb_type_atomlist 1211 use xtb_setparam 1212 implicit none 1213 character(len=*), parameter :: source = 'userdata_reactor' 1214 type(TEnvironment), intent(inout) :: env 1215 character(len=*),intent(in) :: key 1216 character(len=*),intent(in) :: val 1217 integer, intent(in) :: nat 1218 integer, intent(in) :: at(nat) 1219 type(TIdentityMap), intent(in) :: idMap 1220 real(wp),intent(in) :: xyz(3,nat) 1221 1222 type(TAtomList) :: atl 1223 1224 integer :: idum 1225 real(wp) :: ddum 1226 logical :: ldum 1227 integer :: nlist 1228 integer, allocatable :: list(:) 1229 integer :: i,j 1230 1231 integer :: narg 1232 character(len=p_str_length),dimension(p_arg_length) :: argv 1233 1234 call parse(val,comma,argv,narg) 1235 if (verbose) then 1236 do idum = 1, narg 1237 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 1238 enddo 1239 endif 1240 select case(key) 1241 case default ! ignore, don't even think about raising them 1242 case('atoms') 1243 call atl%new(val) 1244 if (atl%get_error()) then 1245 call env%warning('something is wrong in the reactor atom list',source) 1246 return 1247 endif 1248 if (reactset%nat > 0) call atl%add(reactset%atoms(:reactset%nat)) 1249 call atl%to_list(list) 1250 reactset%atoms = list 1251 reactset%nat = size(list) 1252 end select 1253 1254end subroutine set_reactor 1255 1256subroutine set_path(env,key,val,nat,at,idMap,xyz) 1257 use xtb_type_atomlist 1258 use xtb_setparam 1259 implicit none 1260 character(len=*), parameter :: source = 'userdata_path' 1261 type(TEnvironment), intent(inout) :: env 1262 character(len=*),intent(in) :: key 1263 character(len=*),intent(in) :: val 1264 integer, intent(in) :: nat 1265 integer, intent(in) :: at(nat) 1266 type(TIdentityMap), intent(in) :: idMap 1267 real(wp),intent(in) :: xyz(3,nat) 1268 1269 type(TAtomList) :: atl 1270 1271 integer :: idum 1272 real(wp) :: ddum 1273 logical :: ldum 1274 integer :: nlist 1275 integer, allocatable :: list(:) 1276 integer :: i,j 1277 1278 integer :: narg 1279 character(len=p_str_length),dimension(p_arg_length) :: argv 1280 1281 call parse(val,comma,argv,narg) 1282 if (verbose) then 1283 do idum = 1, narg 1284 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 1285 enddo 1286 endif 1287 select case(key) 1288 case default ! ignore, don't even think about raising them 1289 case('atoms') 1290 call atl%new(val) 1291 if (atl%get_error()) then 1292 call env%warning('something is wrong in the bias path atom list',source) 1293 return 1294 endif 1295 if (pathset%nat > 0) call atl%add(pathset%atoms(:pathset%nat)) 1296 call atl%to_list(list) 1297 pathset%atoms = list 1298 pathset%nat = size(list) 1299 end select 1300 1301end subroutine set_path 1302 1303subroutine set_metadyn(env,key,val,nat,at,idMap,xyz) 1304 use xtb_type_atomlist 1305 use xtb_fixparam 1306 implicit none 1307 character(len=*), parameter :: source = 'userdata_metadyn' 1308 type(TEnvironment), intent(inout) :: env 1309 character(len=*),intent(in) :: key 1310 character(len=*),intent(in) :: val 1311 integer, intent(in) :: nat 1312 integer, intent(in) :: at(nat) 1313 type(TIdentityMap), intent(in) :: idMap 1314 real(wp),intent(in) :: xyz(3,nat) 1315 1316 type(TAtomList) :: atl 1317 1318 integer :: idum 1319 real(wp) :: ddum 1320 logical :: ldum 1321 integer :: nlist 1322 integer, allocatable :: list(:) 1323 integer :: i,j,iat 1324 1325 integer :: narg 1326 character(len=p_str_length),dimension(p_arg_length) :: argv 1327 1328 call parse(val,comma,argv,narg) 1329 if (verbose) then 1330 do idum = 1, narg 1331 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 1332 enddo 1333 endif 1334 select case(key) 1335 case default ! ignore, don't even think about raising them 1336 case('bias atoms') 1337 call atl%new(val) 1338 if (atl%get_error()) then 1339 call env%warning('something is wrong in the metadynamic atom list',source) 1340 return 1341 endif 1342 if (rmsdset%nat > 0) call atl%add(rmsdset%atoms(:rmsdset%nat)) 1343 call atl%to_list(list) 1344 rmsdset%atoms = list 1345 rmsdset%nat = size(list) 1346 case('bias elements') 1347 call atl%new 1348 do idum = 1, narg 1349 ! get element by symbol 1350 if (idMap%has(argv(idum))) then 1351 call idMap%get(list, argv(idum)) 1352 if (allocated(list)) then 1353 call atl%add(list) 1354 else 1355 call env%warning("Unknown element: '"//trim(argv(idum))//"'",source) 1356 cycle 1357 end if 1358 else 1359 ldum = getValue(env,trim(argv(idum)),iat) 1360 if (.not.ldum) cycle ! skip garbage input 1361 ! check for unreasonable input 1362 if (iat > 0) then 1363 call atl%add(at.eq.iat) 1364 else 1365 call env%warning("Unknown element: '"//trim(argv(idum))//"'",source) 1366 cycle 1367 endif 1368 endif 1369 enddo 1370 if (rmsdset%nat > 0) call atl%add(rmsdset%atoms(:rmsdset%nat)) 1371 call atl%to_list(list) 1372 rmsdset%atoms = list 1373 rmsdset%nat = size(list) 1374 case('atoms') 1375 call atl%new(val) 1376 if (atl%get_error()) then 1377 call env%warning('something is wrong in the metadynamic atom list',source) 1378 return 1379 endif 1380 if (metaset%nat > 0) call atl%add(metaset%atoms(:metaset%nat)) 1381 call atl%to_list(list) 1382 metaset%atoms = list 1383 metaset%nat = size(list) 1384 case('modify factor') 1385 if (mod(narg,2).ne.0) then 1386 call env%warning("Something went wrong in set_metadyn_ 'modify'.",source) 1387 endif 1388 do i = 1, narg, 2 1389 j = i+1 1390 if (getValue(env,trim(argv(i)),idum).and.& 1391 getValue(env,trim(argv(j)),ddum)) then 1392 if (idum.gt.metaset%maxsave) then 1393 call env%warning('Attempted using factor not present in system.',source) 1394 cycle 1395 endif 1396 metaset%factor(idum) = ddum 1397 endif 1398 enddo 1399 case('scale factor') 1400 if (mod(narg,2).ne.0) then 1401 call env%warning("Something went wrong in set_metadyn_ 'scale'.",source) 1402 endif 1403 do i = 1, narg, 2 1404 j = i+1 1405 if (getValue(env,trim(argv(i)),idum).and.& 1406 getValue(env,trim(argv(j)),ddum)) then 1407 if (idum.gt.metaset%maxsave) then 1408 call env%warning('Attempted scaling factor not present in system.',source) 1409 cycle 1410 endif 1411 metaset%factor(idum) = metaset%factor(idum)*ddum 1412 endif 1413 enddo 1414 end select 1415 1416end subroutine set_metadyn 1417 1418subroutine set_freeze(env,key,val,nat,at,idMap,xyz) 1419 implicit none 1420 character(len=*), parameter :: source = 'userdata_freeze' 1421 type(TEnvironment), intent(inout) :: env 1422 character(len=*),intent(in) :: key 1423 character(len=*),intent(in) :: val 1424 integer, intent(in) :: nat 1425 integer, intent(in) :: at(nat) 1426 type(TIdentityMap), intent(in) :: idMap 1427 real(wp),intent(in) :: xyz(3,nat) 1428 1429 integer :: idum 1430 real(wp) :: ddum 1431 logical :: ldum 1432 integer :: list(nat),nlist 1433 integer :: i,j 1434 1435 integer :: narg 1436 character(len=p_str_length),dimension(p_arg_length) :: argv 1437 1438 call parse(val,comma,argv,narg) 1439 if (verbose) then 1440 do idum = 1, narg 1441 write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum)) 1442 enddo 1443 endif 1444 select case(key) 1445 case default ! ignore, don't even think about raising them 1446 case('hessf') 1447 case('hessa') 1448 end select 1449 1450end subroutine set_freeze 1451 1452subroutine set_legacy(env,key,val,nat,at,idMap,xyz) 1453 implicit none 1454 character(len=*), parameter :: source = 'userdata_legacy' 1455 type(TEnvironment), intent(inout) :: env 1456 character(len=*),intent(in) :: key 1457 character(len=*),intent(in) :: val 1458 integer, intent(in) :: nat 1459 integer, intent(in) :: at(nat) 1460 type(TIdentityMap), intent(in) :: idMap 1461 real(wp),intent(in) :: xyz(3,nat) 1462 integer :: err 1463 integer :: idum 1464 real(wp) :: ddum 1465 logical :: ldum 1466 select case(key) 1467 case default ! complaining about unknown keywords should already be done 1468 continue ! so we do nothing here 1469 case('hessf'); call set_fix(env,'freeze',val,nat,at,idMap,xyz) 1470! case('hessa'); call set_frozh('hessa',val) 1471 case('fragment1'); call set_split(env,'fragment1',val,nat,at,idMap,xyz) 1472 case('fragment2'); call set_split(env,'fragment1',val,nat,at,idMap,xyz) 1473 case('constrxyz'); call set_fix(env,'atoms',val,nat,at,idMap,xyz) 1474! case('constrainel') 1475! case('constrain') 1476! case('scan') 1477 case('ellips'); call set_wall(env,'ellipsoid',val,nat,at,idMap,xyz) 1478 case('sphere'); call set_wall(env,'sphere',val,nat,at,idMap,xyz) 1479 case('fix'); call set_fix(env,'atoms',val,nat,at,idMap,xyz) 1480 case('atomlist+'); call set_metadyn(env,'atoms',val,nat,at,idMap,xyz) 1481 end select 1482 1483end subroutine set_legacy 1484 1485end module xtb_constrain_param 1486