1c 2c $Id$ 3c 4c 5c This file contains a list of beads, where a bead is defined 6c as a movecs_filename, geom_name, and a Energy. This list is used 7c in the NEB code. 8c 9c Author - Eric Bylaska 10c Creation Date - 11/11/2002 11c 12c 13c *************************************************** 14c * * 15c * create_bead_list * 16c * * 17c *************************************************** 18 19c This routine initializes a bead_list labeled 'tag' on 20c the runtime database, rtdb. 21c 22c Entry - rtdb: pointeer to runtime database 23c tag: name of created bead list 24c perm_movecs_name: name of movecs files used by call to task 25c 26c 27 subroutine init_bead_list(rtdb,tag,perm_movecs_name) 28 implicit none 29 integer rtdb 30 character*(*) tag,perm_movecs_name 31 32#include "rtdb.fh" 33#include "mafdecls.fh" 34 35* ***** local variables **** 36 logical value 37 integer size,taglen,movecslen 38 character*255 rtdb_name 39 40 integer bead_rtdb 41 common /bead_list/ bead_rtdb 42 43* **** external functions **** 44 integer inp_strlen 45 external inp_strlen 46 47 bead_rtdb = rtdb 48 49 taglen = inp_strlen(tag) 50 movecslen = inp_strlen(perm_movecs_name) 51 52 size = 0 53 rtdb_name = tag(1:taglen)//':size' 54 value = rtdb_put(bead_rtdb,rtdb_name,mt_int,1,size) 55 56 rtdb_name = tag(1:taglen)//':perm_movecs' 57 value = value.and. 58 > rtdb_cput(bead_rtdb,rtdb_name,1, 59 > perm_movecs_name(1:movecslen)) 60 61 if (.not.value) call errquit('init_bead_list failed',0,0) 62 63 return 64 end 65 66 subroutine set_rtdb_bead_list(rtdb) 67 implicit none 68 integer rtdb 69 70 integer bead_rtdb 71 common /bead_list/ bead_rtdb 72 73 bead_rtdb = rtdb 74 return 75 end 76 77 78 subroutine reset_bead_list(tag) 79 implicit none 80 integer rtdb 81 character*(*) tag 82 83#include "rtdb.fh" 84#include "mafdecls.fh" 85 86* ***** local variables **** 87 logical value 88 integer size,taglen 89 character*255 rtdb_name 90 91 integer bead_rtdb 92 common /bead_list/ bead_rtdb 93 94* **** external functions **** 95 integer inp_strlen 96 external inp_strlen 97 98 taglen = inp_strlen(tag) 99 size = 0 100 rtdb_name = tag(1:taglen)//':size' 101 value = rtdb_put(bead_rtdb,rtdb_name,mt_int,1,size) 102 103 return 104 end 105 106c *************************************************** 107c * * 108c * add_bead_list * 109c * * 110c *************************************************** 111 112c This routine add a bead to the bead_list labeled 'tag'. 113c 114c Entry - tag: name of created bead list 115c movecs_name: movecs filename of added bead 116c geom_name: geometry name of added bead 117c 118c 119 subroutine add_bead_list(tag,movecs_name,geom_name) 120 implicit none 121 character*(*) tag,movecs_name,geom_name 122 123#include "rtdb.fh" 124#include "mafdecls.fh" 125 126* ***** local variables **** 127 logical value 128 integer size,taglen,movecslen,geomlen 129 character*255 rtdb_name 130 real*8 stress(9) 131 132 integer bead_rtdb 133 common /bead_list/ bead_rtdb 134 135* **** external functions **** 136 character*7 bead_index_name 137 integer inp_strlen 138 external bead_index_name 139 external inp_strlen 140 141 142 taglen = inp_strlen(tag) 143 movecslen = inp_strlen(movecs_name) 144 geomlen = inp_strlen(geom_name) 145 146 147 rtdb_name = tag(1:taglen)//':size' 148 value = rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size) 149 size = size+1 150 value = value.and. 151 > rtdb_put(bead_rtdb,rtdb_name,mt_int,1,size) 152 153 154 rtdb_name = tag(1:taglen)//bead_index_name(size)//':movecs_name' 155 value = value.and. 156 > rtdb_cput(bead_rtdb,rtdb_name,1,movecs_name(1:movecslen)) 157 158 rtdb_name = tag(1:taglen)//bead_index_name(size)//':geom_name' 159 value = value.and. 160 > rtdb_cput(bead_rtdb,rtdb_name,1,geom_name(1:geomlen)) 161 162 rtdb_name = tag(1:taglen)//bead_index_name(size)//':energy' 163 value = value.and. 164 > rtdb_put(bead_rtdb,rtdb_name,mt_dbl,1,0.0d0) 165 166 call dcopy(9,0.0d0,0,stress,1) 167 rtdb_name = tag(1:taglen)//bead_index_name(size)//':stress' 168 value = value.and. 169 > rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,stress) 170 171 if (.not.value) call errquit('add_bead_list failed',0,0) 172 173 return 174 end 175 176 177 178c *************************************************** 179c * * 180c * delete_bead_list * 181c * * 182c *************************************************** 183 184c This routine deletes bead,i, from the bead_list labeled 'tag'. 185c 186c Entry - tag: name of created bead list 187c movecs_name: movecs filename of added bead 188c geom_name: geometry name of added bead 189c 190c 191 subroutine delete_bead_list(tag,i) 192 implicit none 193 character*(*) tag 194 integer i 195 196#include "rtdb.fh" 197#include "mafdecls.fh" 198 199* ***** local variables **** 200 logical value 201 integer size,taglen,tmplen,j 202 real*8 energy,stress(9) 203 character*255 rtdb_name,tmp_name 204 205 integer bead_rtdb 206 common /bead_list/ bead_rtdb 207 208* **** external functions **** 209 character*7 bead_index_name 210 integer inp_strlen 211 external bead_index_name 212 external inp_strlen 213 214 215 taglen = inp_strlen(tag) 216 217 rtdb_name = tag(1:taglen)//':size' 218 value = rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size) 219 if (i.gt.size) return 220 221 value = value.and.rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size-1) 222 223 do j=i,size-1 224 225 rtdb_name = tag(1:taglen)//bead_index_name(j+1)//':movecs_name' 226c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name)) 227c > tmp_name = 'bead'//bead_index_name(j+1)//'.movecs' 228 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name)) 229 > tmp_name = tag(1:taglen)//bead_index_name(j+1)//'.movecs' 230 tmplen = inp_strlen(tmp_name) 231 rtdb_name = tag(1:taglen)//bead_index_name(j)//':movecs_name' 232 value = value.and. 233 > rtdb_cput(bead_rtdb,rtdb_name,1,tmp_name(1:tmplen)) 234 235 rtdb_name = tag(1:taglen)//bead_index_name(j+1)//':geom_name' 236c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name)) 237c > tmp_name = 'bead'//bead_index_name(j+1)//':geom' 238 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name)) 239 > tmp_name = tag(1:taglen)//bead_index_name(j+1)//':geom' 240 tmplen = inp_strlen(tmp_name) 241 rtdb_name = tag(1:taglen)//bead_index_name(j)//':geom_name' 242 value = value.and. 243 > rtdb_cput(bead_rtdb,rtdb_name,1,tmp_name(1:tmplen)) 244 245 rtdb_name = tag(1:taglen)//bead_index_name(j+1)//':energy' 246 value = value.and. 247 > rtdb_get(bead_rtdb,rtdb_name,mt_dbl,1,energy) 248 rtdb_name = tag(1:taglen)//bead_index_name(j)//':energy' 249 value = value.and. 250 > rtdb_put(bead_rtdb,rtdb_name,mt_dbl,1,energy) 251 252 rtdb_name = tag(1:taglen)//bead_index_name(j+1)//':stress' 253 value = value.and. 254 > rtdb_get(bead_rtdb,rtdb_name,mt_dbl,9,stress) 255 rtdb_name = tag(1:taglen)//bead_index_name(j)//':stress' 256 value = value.and. 257 > rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,stress) 258 259 end do 260 261 if (.not.value) call errquit('delete_bead_list failed',0,0) 262 263 return 264 end 265 266 267c *************************************************** 268c * * 269c * size_bead_list * 270c * * 271c *************************************************** 272 273c This routine returns the number of bead in the bead_list 274c labeled "tag". 275c 276c Entry - tag: name of bead list 277c Exit - returns number of beads 278c 279 integer function size_bead_list(tag) 280 implicit none 281 character*(*) tag 282 283#include "rtdb.fh" 284#include "mafdecls.fh" 285 286* ***** local variables **** 287 logical value 288 integer size,taglen 289 character*255 rtdb_name 290 291 integer bead_rtdb 292 common /bead_list/ bead_rtdb 293 294* **** external functions **** 295 integer inp_strlen 296 external inp_strlen 297 298 taglen = inp_strlen(tag) 299 300 rtdb_name = tag(1:taglen)//':size' 301 value = rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size) 302 if (.not.value) call errquit('size_bead_list failed',0,0) 303 304 size_bead_list = size 305 return 306 end 307 308 309c *************************************************** 310c * * 311c * nion_bead_list * 312c * * 313c *************************************************** 314 315c This routine returns the number of ions, nions, in bead,i, 316c from the bead_list labeled 'tag'. 317c 318c Entry - tag: name of bead list 319c i: bead number 320c Exit - returns the number ions in bead, i 321c 322 integer function nion_bead_list(tag,i) 323 implicit none 324 character*(*) tag 325 integer i 326 327#include "rtdb.fh" 328#include "mafdecls.fh" 329#include "geom.fh" 330 331* ***** local variables **** 332 logical value,ostress 333 integer size,taglen,nion,geom,geomlen 334 character*255 rtdb_name,geom_name 335 336 integer bead_rtdb 337 common /bead_list/ bead_rtdb 338 339* **** external functions **** 340 logical bead_includestress 341 integer inp_strlen,size_bead_list 342 character*7 bead_index_name 343 external bead_includestress 344 external inp_strlen,size_bead_list 345 external bead_index_name 346 347 taglen = inp_strlen(tag) 348 349 size = size_bead_list(tag) 350 if (i.gt.size) then 351 nion_bead_list = 0 352 return 353 end if 354 355 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 356c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 357c > geom_name = 'bead'//bead_index_name(i)//':geom' 358 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 359 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 360 geomlen = inp_strlen(geom_name) 361 362 value = geom_create(geom,'bead_tmp') 363 value = value.and.geom_rtdb_load(bead_rtdb,geom, 364 > geom_name(1:geomlen)) 365 value = value.and.geom_ncent(geom,nion) 366 value = value.and.geom_destroy(geom) 367 if (.not.value) call errquit('nion_bead_list failed',0,0) 368 369 if (bead_includestress()) nion = nion + 3 370 371 nion_bead_list = nion 372 return 373 end 374 375 376c *************************************************** 377c * * 378c * run_bead_list * 379c * * 380c *************************************************** 381 382c This routine runs bead,i, from the bead_list labeled "tag". 383c 384c Entry - tag: name of bead list 385c i: bead number 386c 387c 388 subroutine run_bead_list(tag,i,task_gradient_gen) 389 implicit none 390 character*(*) tag 391 integer i 392 393#include "rtdb.fh" 394#include "mafdecls.fh" 395#include "geom.fh" 396#include "global.fh" 397#include "errquit.fh" 398 399* ***** local variables **** 400 logical value,status,ignore 401 integer geom,gradient(2),nion,nelem 402 integer size,taglen,permlen,movecslen,geomlen 403 real*8 energy,stress(9) 404 character*255 rtdb_name,perm_name,movecs_name,geom_name 405 character*32 theory 406 407 integer bead_rtdb 408 common /bead_list/ bead_rtdb 409 410* **** external functions **** 411 logical task_gradient_gen,bead_includestress 412 external task_gradient_gen,bead_includestress 413 character*7 bead_index_name 414 integer inp_strlen 415 external bead_index_name 416 external inp_strlen 417 418 value = .true. 419 taglen = inp_strlen(tag) 420 421 rtdb_name = tag(1:taglen)//':perm_movecs' 422 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,perm_name)) 423 > call util_file_prefix('movecs',perm_name) 424 call util_file_name_resolve(perm_name, .false.) 425 permlen = inp_strlen(perm_name) 426 427 rtdb_name = tag(1:taglen)//':size' 428 if (.not.rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size)) 429 > size = 0 430 if (i.gt.size) return 431 432 rtdb_name = tag(1:taglen)//bead_index_name(i)//':movecs_name' 433c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,movecs_name)) 434c > movecs_name = 'bead'//bead_index_name(i)//'.movecs' 435 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,movecs_name)) 436 > movecs_name = tag(1:taglen)//bead_index_name(i)//'.movecs' 437 438 call util_file_name_resolve(movecs_name, .false.) 439 movecslen = inp_strlen(movecs_name) 440 441 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 442c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 443c > geom_name = 'bead'//bead_index_name(i)//':geom' 444 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 445 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 446 geomlen = inp_strlen(geom_name) 447 448 value = value.and. 449 > rtdb_cput(bead_rtdb,'geometry',1,geom_name(1:geomlen)) 450 451 if (.not.value) call errquit('run_bead_list failed',0,0) 452 453 454* *** copy bead movecs to taskmovecs *** 455 if (ga_nodeid().eq.0) then 456 inquire(file=movecs_name,exist=status) 457 if (status) then 458 call util_file_copy(movecs_name(1:movecslen), 459 > perm_name(1:permlen)) 460 end if 461 end if 462 463* *** run gradient task *** 464 if (.not.rtdb_get(bead_rtdb,"task:ignore",mt_log,1,ignore)) 465 > ignore = .false. 466 467 if(.not.rtdb_put(bead_rtdb,"scf:converged",mt_log,1,.false.)) 468 > call errquit('scf:converged put',0,0) 469 if(.not.rtdb_put(bead_rtdb,"dft:converged",mt_log,1,.false.)) 470 > call errquit('dft:converged put',0,0) 471 472 if(.not.rtdb_put(bead_rtdb,"neb:ibead",mt_int,1,i)) 473 > call errquit('neb:ibead put',0,0) 474 if(ga_nodeid().eq.0) write(*,*) "neb: running bead ",i 475 value=value.and.(task_gradient_gen(bead_rtdb).or.ignore) 476 if ((.not.value)) 477 > call errquit('run_bead_list failed',1,0) 478 479 480* *** copy taskmovecs to bead movecs **** 481 if (ga_nodeid().eq.0) then 482 inquire(file=perm_name,exist=status) 483 if (status) then 484 call util_file_copy(perm_name(1:permlen), 485 > movecs_name(1:movecslen)) 486 end if 487 end if 488 489* ***** set the energy ***** 490 value = value.and. 491 > rtdb_get(bead_rtdb,'task:energy',mt_dbl,1,energy) 492 rtdb_name = tag(1:taglen)//bead_index_name(i)//':energy' 493 value = value.and. 494 > rtdb_put(bead_rtdb,rtdb_name,mt_dbl,1,energy) 495 if (.not.value) call errquit('run_bead_list failed',2,0) 496 497 498* ***** set the forces ****** 499 value = value.and.geom_create(geom,'geometry') 500 value = value.and.geom_rtdb_load(bead_rtdb,geom,'geometry') 501 value = value.and.geom_ncent(geom,nion) 502 if (.not.value) call errquit('run_bead_list failed',3,0) 503 504 value = value.and. 505 > MA_push_get(mt_dbl,(3*nion), 506 > 'gradient',gradient(2),gradient(1)) 507 value = value.and. 508 > rtdb_get(bead_rtdb,'task:gradient',mt_dbl,(3*nion), 509 > dbl_mb(gradient(1))) 510 511 value = value.and.geom_vel_set(geom,dbl_mb(gradient(1))) 512 value = value.and.MA_pop_stack(gradient(2)) 513 value = value.and.geom_rtdb_delete(bead_rtdb,'geometry') 514 value = value.and.geom_rtdb_store(bead_rtdb,geom,'geometry') 515 value = value.and.geom_destroy(geom) 516 if (.not.value) call errquit('run_bead_list failed',4,0) 517 518 519* ***** set the stresses ****** 520 if (bead_includestress()) then 521 522 if (.not.rtdb_cget(bead_rtdb, 'task:theory', 1, theory)) 523 > call errquit('run_bead_list: stress theory not specified', 524 > 0,RTDB_ERR) 525 if (theory.eq.'pspw') then 526 if (.not.rtdb_get(bead_rtdb, 'pspw:stress', mt_dbl, 9,stress)) 527 > call errquit('run_bead_list: could not get stress',0,0) 528 else if (theory.eq.'band') then 529 if (.not.rtdb_get(bead_rtdb, 'band:stress', mt_dbl, 9,stress)) 530 > call errquit('run_bead_list: could not get stress',0,0) 531 else if (theory.eq.'paw') then 532 if (.not.rtdb_get(bead_rtdb, 'paw:stress', mt_dbl, 9,stress)) 533 > call errquit('run_bead_list: could not get stress',0,0) 534 else 535 call errquit('run_bead_list: no stress in theory',0,RTDB_ERR) 536 end if 537 538 rtdb_name = tag(1:taglen)//bead_index_name(i)//':stress' 539 value = value.and. 540 > rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,stress) 541 if (.not.value) call errquit('run_bead_list failed',2,0) 542 end if 543 544 if(ga_nodeid().eq.0) then 545 write(*,*) "neb: finished bead ",i 546 write(*,*) "neb: final energy",energy 547 end if 548 549 return 550 end 551 552 553c *************************************************** 554c * * 555c * runall_bead_list * 556c * * 557c *************************************************** 558 559c This routine runs all beads from the bead_list labeled 'tag'. 560c 561c Entry - tag: name of bead list 562c 563c 564 subroutine runall_bead_list(tag,task_gradient_gen) 565 implicit none 566 character*(*) tag 567 logical task_gradient_gen 568 external task_gradient_gen 569 570* *** local variables *** 571 integer i,nbeads 572 573* *** external functions *** 574 integer size_bead_list 575 external size_bead_list 576 577 nbeads = size_bead_list(tag) 578 do i=1,nbeads 579 call run_bead_list(tag,i,task_gradient_gen) 580 end do 581 582 return 583 end 584 585 586 587c *************************************************** 588c * * 589c * runmid_bead_list * 590c * * 591c *************************************************** 592 593c This routine runs all beads between 2 and nbeads-1 594c from the bead_list labeled 'tag'. 595c 596c Entry - tag: name of bead list 597c 598c 599 subroutine runmid_bead_list(tag,task_gradient_gen, 600 > freeze1,freezen) 601 implicit none 602 character*(*) tag 603 logical task_gradient_gen 604 external task_gradient_gen 605 logical freeze1,freezen 606 607* *** local variables *** 608 integer i,nbeads,istart,iend 609 610* *** external functions *** 611 integer size_bead_list 612 external size_bead_list 613 614 nbeads = size_bead_list(tag) 615 !do i=2,(nbeads-1) 616 istart = 2 617 iend = nbeads-1 618 if (.not.freeze1) istart = 1 619 if (.not.freezen) iend = nbeads 620 do i=istart,iend 621 call run_bead_list(tag,i,task_gradient_gen) 622 end do 623 624 return 625 end 626 627 628 629c *************************************************** 630c * * 631c * energy_bead_list * 632c * * 633c *************************************************** 634 635c This routine returns the current energy of bead, i, 636c from the bead_list labeled 'tag'. 637c 638c Entry - tag: name of bead list 639c i: bead number 640c 641 real*8 function energy_bead_list(tag,i) 642 implicit none 643 character*(*) tag 644 integer i 645 646#include "rtdb.fh" 647#include "mafdecls.fh" 648 649* ***** local variables **** 650 logical value 651 integer taglen 652 real*8 energy 653 character*255 rtdb_name 654 655 integer bead_rtdb 656 common /bead_list/ bead_rtdb 657 658* **** external functions **** 659 integer inp_strlen 660 character*7 bead_index_name 661 external inp_strlen 662 external bead_index_name 663 664 taglen = inp_strlen(tag) 665 666 rtdb_name = tag(1:taglen)//bead_index_name(i)//':energy' 667 value = rtdb_get(bead_rtdb,rtdb_name,mt_dbl,1,energy) 668 if (.not.value) call errquit('energy_bead_list failed',0,0) 669 670 energy_bead_list = energy 671 return 672 end 673 674 675 676c *************************************************** 677c * * 678c * coords_get_bead_list * 679c * * 680c *************************************************** 681 682c This routine returns the ion coordinates, c, of the "i" bead 683c from the "tag" bead list. 684c 685c Entry - tag: name of bead list 686c i: bead number 687c Exit - c: ion coordinates 688c 689 subroutine coords_get_bead_list(tag,i,c) 690 implicit none 691 character*(*) tag 692 integer i 693 real*8 c(3,*) 694 695#include "rtdb.fh" 696#include "mafdecls.fh" 697#include "geom.fh" 698 699* ***** local variables **** 700 logical value 701 integer geom,nion 702 integer size,taglen,geomlen 703 real*8 energy 704 character*255 rtdb_name,geom_name 705 706 integer bead_rtdb 707 common /bead_list/ bead_rtdb 708 709* **** external functions **** 710 logical bead_includestress,geom_grad_cart_to_frac 711 character*7 bead_index_name 712 integer inp_strlen,size_bead_list 713 external bead_includestress,geom_grad_cart_to_frac 714 external bead_index_name 715 external inp_strlen,size_bead_list 716 717 size = size_bead_list(tag) 718 if (i.gt.size) return 719 720 taglen = inp_strlen(tag) 721 722 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 723c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 724c > geom_name = 'bead'//bead_index_name(i)//':geom' 725 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 726 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 727 geomlen = inp_strlen(geom_name) 728 729 value = geom_create(geom,'bead_tmp') 730 value = value.and.geom_rtdb_load(bead_rtdb,geom, 731 > geom_name(1:geomlen)) 732 value = value.and.geom_ncent(geom,nion) 733 value = value.and.geom_cart_coords_get(geom,c) 734 735 if (bead_includestress()) then 736 737* **** put gradient into fractional **** 738 if (.not. geom_grad_cart_to_frac(geom,c)) 739 $ call errquit('coords_get_bead_list: cart_to_frac?',0,0) 740 741* **** get stress part of gradient *** 742 if (.not. geom_amatrix_get(geom,c(1,nion+1))) 743 > call errquit('coords_get_bead_list:failed to get amatrix',0,0) 744 end if 745 746 value = value.and.geom_destroy(geom) 747 if (.not.value) call errquit('coords_get_bead_list failed',0,0) 748 749 return 750 end 751 752c *************************************************** 753c * * 754c * coords_set_bead_list * 755c * * 756c *************************************************** 757 758c This routine sets the ion coordinates, c, of "i" bead 759c from the "tag" bead list. 760c 761c Entry - tag: name of bead list 762c i: bead number 763c c: ion coordinates 764c 765c 766 subroutine coords_set_bead_list(tag,i,c) 767 implicit none 768 character*(*) tag 769 integer i 770 real*8 c(3,*) 771 772#include "rtdb.fh" 773#include "mafdecls.fh" 774#include "geom.fh" 775 776* ***** local variables **** 777 logical value 778 integer geom,nion 779 integer size,taglen,geomlen 780 real*8 energy 781 character*255 rtdb_name,geom_name 782 783 integer bead_rtdb 784 common /bead_list/ bead_rtdb 785 786* **** external functions **** 787 logical bead_includestress,geom_amatrix_set 788 character*7 bead_index_name 789 integer inp_strlen,size_bead_list 790 external bead_includestress,geom_amatrix_set 791 external bead_index_name 792 external inp_strlen,size_bead_list 793 794 size = size_bead_list(tag) 795 if (i.gt.size) return 796 797 taglen = inp_strlen(tag) 798 799 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 800c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 801c > geom_name = 'bead'//bead_index_name(i)//':geom' 802 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 803 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 804 geomlen = inp_strlen(geom_name) 805 806 value = geom_create(geom,'bead_tmp') 807 value = value.and.geom_rtdb_load(bead_rtdb,geom, 808 > geom_name(1:geomlen)) 809 value = value.and.geom_ncent(geom,nion) 810 811 if (bead_includestress()) then 812 813 if (.not. geom_frac_to_cart(geom, c)) 814 $ call errquit('coords_set_bead_list: frac_to_cart?',0,0) 815 816 if (.not. geom_amatrix_set(geom,c(1,nion+1))) 817 $ call errquit('coords_set_bead_list:failed to set amatrix',0,0) 818 819 820 end if 821 822 value = value.and.geom_cart_coords_set(geom,c) 823 value = value.and.geom_rtdb_delete(bead_rtdb,geom_name(1:geomlen)) 824 value = value.and.geom_rtdb_store(bead_rtdb,geom, 825 > geom_name(1:geomlen)) 826 value = value.and.geom_destroy(geom) 827 if (.not.value) call errquit('coords_set_bead_list failed',0,0) 828 829 return 830 end 831 832 833 834c ********************************************** 835c * * 836c * systype_bead_list * 837c * * 838c ********************************************** 839 840c This routine returns if the systype, i.e. the number of dimensions 841c from the "tag" bead list. Uses bead index=1. 842c 843c Entry - tag: name of bead list 844c Exit - 845c 846 integer function systype_bead_list(tag) 847 implicit none 848 character*(*) tag 849 850#include "rtdb.fh" 851#include "mafdecls.fh" 852#include "geom.fh" 853 854* ***** local variables **** 855 integer i,isystype 856 real*8 amat(3,3) 857 logical value 858 integer geom,nion 859 integer size,taglen,geomlen 860 real*8 energy 861 character*255 rtdb_name,geom_name 862 863 integer bead_rtdb 864 common /bead_list/ bead_rtdb 865 866* **** external functions **** 867 character*7 bead_index_name 868 integer inp_strlen,size_bead_list 869 external bead_index_name 870 external inp_strlen,size_bead_list 871 872 isystype = 0 873 i = 1 874 size = size_bead_list(tag) 875 if (i.gt.size) then 876 systype_bead_list = isystype 877 return 878 end if 879 880 taglen = inp_strlen(tag) 881 882 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 883c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 884c > geom_name = 'bead'//bead_index_name(i)//':geom' 885 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 886 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 887 geomlen = inp_strlen(geom_name) 888 889 value = geom_create(geom,'bead_tmp') 890 value = value.and.geom_rtdb_load(bead_rtdb,geom, 891 > geom_name(1:geomlen)) 892 893* **** get systype *** 894 if (.not.geom_systype_get(geom,isystype)) then 895 isystype = 0 896 end if 897 898 value = value.and.geom_destroy(geom) 899 if (.not.value) call errquit('systype_bead_list failed',0,0) 900 901 systype_bead_list = isystype 902 return 903 end 904 905 906 907 908c *************************************************** 909c * * 910c * amatrix_get_bead_list * 911c * * 912c *************************************************** 913 914c This routine returns the amatrix, amat, of the "i" bead 915c from the "tag" bead list. 916c 917c Entry - tag: name of bead list 918c i: bead number 919c Exit - amat: amatrix 920c 921 subroutine amatrix_get_bead_list(tag,i,amat) 922 implicit none 923 character*(*) tag 924 integer i 925 real*8 amat(3,3) 926 927#include "rtdb.fh" 928#include "mafdecls.fh" 929#include "geom.fh" 930 931* ***** local variables **** 932 logical value 933 integer geom,nion 934 integer size,taglen,geomlen 935 real*8 energy 936 character*255 rtdb_name,geom_name 937 938 integer bead_rtdb 939 common /bead_list/ bead_rtdb 940 941* **** external functions **** 942 character*7 bead_index_name 943 integer inp_strlen,size_bead_list 944 external bead_index_name 945 external inp_strlen,size_bead_list 946 947 size = size_bead_list(tag) 948 if (i.gt.size) return 949 950 taglen = inp_strlen(tag) 951 952 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 953c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 954c > geom_name = 'bead'//bead_index_name(i)//':geom' 955 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 956 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 957 geomlen = inp_strlen(geom_name) 958 959 value = geom_create(geom,'bead_tmp') 960 value = value.and.geom_rtdb_load(bead_rtdb,geom, 961 > geom_name(1:geomlen)) 962 963* **** get amatrix *** 964 if (.not.geom_amatrix_get(geom,amat)) then 965 call dcopy(9,0.0d0,0,amat,1) 966 end if 967 968 value = value.and.geom_destroy(geom) 969 if (.not.value) call errquit('amatrix_get_bead_list failed',0,0) 970 971 return 972 end 973 974 975 976 977c *************************************************** 978c * * 979c * gradient_get_bead_list * 980c * * 981c *************************************************** 982 983c This routine returns the ion gradients, c, of "i" bead from 984c the "tag" bead list. 985c 986c Entry - tag: name of bead list 987c i: bead number 988c Exit - c: ion forces 989c 990 subroutine gradient_get_bead_list(tag,i,c) 991 implicit none 992 character*(*) tag 993 integer i 994 real*8 c(3,*) 995 996#include "rtdb.fh" 997#include "mafdecls.fh" 998#include "geom.fh" 999 1000* ***** local variables **** 1001 logical value 1002 integer geom,nion 1003 integer size,taglen,geomlen 1004 real*8 energy 1005 character*255 rtdb_name,geom_name 1006 1007 integer bead_rtdb 1008 common /bead_list/ bead_rtdb 1009 1010* **** external functions **** 1011 logical bead_includestress 1012 character*7 bead_index_name 1013 integer inp_strlen,size_bead_list 1014 external bead_includestress 1015 external bead_index_name 1016 external inp_strlen,size_bead_list 1017 1018 size = size_bead_list(tag) 1019 if (i.gt.size) return 1020 1021 taglen = inp_strlen(tag) 1022 1023 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 1024c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 1025c > geom_name = 'bead'//bead_index_name(i)//':geom' 1026 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 1027 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 1028 geomlen = inp_strlen(geom_name) 1029 1030 value = geom_create(geom,'bead_tmp') 1031 value = value.and.geom_rtdb_load(bead_rtdb,geom, 1032 > geom_name(1:geomlen)) 1033 value = value.and.geom_ncent(geom,nion) 1034 value = value.and.geom_vel_get(geom,c) 1035 value = value.and.geom_destroy(geom) 1036 if (.not.value) call errquit('gradient_get_bead_list failed',0,0) 1037 1038 if (bead_includestress()) then 1039 rtdb_name = tag(1:taglen)//bead_index_name(i)//':stress' 1040 value = value.and. 1041 > rtdb_get(bead_rtdb,rtdb_name,mt_dbl,9,c(1,nion+1)) 1042 end if 1043 1044 return 1045 end 1046 1047c *************************************************** 1048c * * 1049c * gradient_set_bead_list * 1050c * * 1051c *************************************************** 1052 1053c This routine sets the ion gradients, c, of "i" bead from the 1054c "tag" bead list'. 1055c 1056c Entry - tag: name of bead list 1057c i: bead number 1058c c: ion forces 1059c 1060c 1061 subroutine gradient_set_bead_list(tag,i,c) 1062 implicit none 1063 character*(*) tag 1064 integer i 1065 real*8 c(3,*) 1066 1067#include "rtdb.fh" 1068#include "mafdecls.fh" 1069#include "geom.fh" 1070 1071* ***** local variables **** 1072 logical value 1073 integer geom,nion 1074 integer size,taglen,geomlen 1075 real*8 energy 1076 character*255 rtdb_name,geom_name 1077 1078 integer bead_rtdb 1079 common /bead_list/ bead_rtdb 1080 1081* **** external functions **** 1082 logical bead_includestress 1083 character*7 bead_index_name 1084 integer inp_strlen,size_bead_list 1085 external bead_includestress 1086 external bead_index_name 1087 external inp_strlen,size_bead_list 1088 1089 size = size_bead_list(tag) 1090 if (i.gt.size) return 1091 1092 taglen = inp_strlen(tag) 1093 1094 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 1095c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 1096c > geom_name = 'bead'//bead_index_name(i)//':geom' 1097 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 1098 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 1099 geomlen = inp_strlen(geom_name) 1100 1101 value = geom_create(geom,'bead_tmp') 1102 value = value.and.geom_rtdb_load(bead_rtdb, 1103 > geom,geom_name(1:geomlen)) 1104 value = value.and.geom_ncent(geom,nion) 1105 value = value.and.geom_vel_set(geom,c) 1106 value = value.and.geom_rtdb_delete(bead_rtdb, 1107 > geom_name(1:geomlen)) 1108 value = value.and.geom_rtdb_store(bead_rtdb, 1109 > geom,geom_name(1:geomlen)) 1110 value = value.and.geom_destroy(geom) 1111 1112 if (bead_includestress()) then 1113 rtdb_name = tag(1:taglen)//bead_index_name(i)//':stress' 1114 value = value.and. 1115 > rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,c(1,nion+1)) 1116 1117 end if 1118 if (.not.value) call errquit('gradient_get_bead_list failed',0,0) 1119 1120 return 1121 end 1122 1123 1124 1125 1126 1127c *************************************************** 1128c * * 1129c * create_xyz_file_bead_list * 1130c * * 1131c *************************************************** 1132 1133c This routine returns the current energy of bead, i, 1134c from the bead_list labeled 'tag'. 1135c 1136c Entry - tag: name of bead list 1137c 1138c 1139 subroutine create_xyz_file_bead_list(luout,tag,header) 1140 implicit none 1141 integer luout 1142 character*(*) tag 1143 logical header 1144 1145#include "rtdb.fh" 1146#include "mafdecls.fh" 1147#include "geom.fh" 1148#include "global.fh" 1149 1150* ***** local variables **** 1151 logical value,isperiodic 1152 integer taglen,i,ii,geom,geomlen,nbeads,nion,i1,i2 1153 real*8 energy,q,rxyz(3),amat(3,3),rr 1154 character*255 rtdb_name,geom_name 1155 character*2 symbol 1156 character*16 t,name 1157 1158 integer bead_rtdb 1159 common /bead_list/ bead_rtdb 1160 1161* **** external functions **** 1162 integer inp_strlen,size_bead_list 1163 character*7 bead_index_name 1164 external inp_strlen,size_bead_list 1165 external bead_index_name 1166 1167 1168 if ((ga_nodeid().eq.0).and.(header)) then 1169 write(luout,*) 1170 write(luout,*) 'XYZ FILE for bead_list:',tag 1171 write(luout,*) '------------------------------------------' 1172 end if 1173 taglen = inp_strlen(tag) 1174 nbeads = size_bead_list(tag) 1175 value = geom_create(geom,'bead_tmp') 1176 do i=1,nbeads 1177 rtdb_name = tag(1:taglen)//bead_index_name(i)//':energy' 1178 value = value.and.rtdb_get(bead_rtdb,rtdb_name, 1179 > mt_dbl,1,energy) 1180 1181 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 1182c if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 1183c > geom_name = 'bead'//bead_index_name(i)//':geom' 1184 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 1185 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 1186 geomlen = inp_strlen(geom_name) 1187 value = value.and.geom_rtdb_load(bead_rtdb, 1188 > geom,geom_name(1:geomlen)) 1189 value = value.and.geom_ncent(geom,nion) 1190 1191 call amatrix_get_bead_list(tag,i,amat) 1192 rr = 0.0d0 1193 do i2=1,3 1194 do i1=1,3 1195 rr = rr + dabs(amat(i1,i2)) 1196 end do 1197 end do 1198 isperiodic = (rr.gt.3.001d0) 1199 1200 1201 if (ga_nodeid().eq.0) then 1202 write(luout,*) nion 1203 if (isperiodic) then 1204 write(luout,'(A,F14.6,A,9F12.6)') 'energy=',energy, 1205 > ' amatrix=',((amat(i1,i2),i1=1,3),i2=1,3) 1206 else 1207 write(luout,'(A,F14.6)') 'energy=',energy 1208 end if 1209 end if 1210 do ii=1,nion 1211 value = value.and.geom_cent_get(geom,ii,t,rxyz,q) 1212 value = value.and. geom_tag_to_element(t,symbol,name,q) 1213 rxyz(1)= rxyz(1)*0.529177d0 1214 rxyz(2)= rxyz(2)*0.529177d0 1215 rxyz(3)= rxyz(3)*0.529177d0 1216 if (ga_nodeid().eq.0) then 1217 write(luout,'(A2,6x,3F12.6)') symbol,rxyz 1218 end if 1219 end do 1220 end do 1221 value = value.and.geom_destroy(geom) 1222 1223 if ((ga_nodeid().eq.0).and.(header)) then 1224 write(luout,*) 1225 end if 1226 1227 if (.not.value) 1228 > call errquit('create_xyz_file_bead_list failed',0,0) 1229 1230 return 1231 end 1232 1233c *************************************************** 1234c * * 1235c * create_cif_file_bead_list * 1236c * * 1237c *************************************************** 1238 1239c This routine returns the current energy and coordinates of bead, i, 1240c from the bead_list labeled 'tag'. 1241c 1242c Entry - tag: name of bead list 1243c 1244c 1245 subroutine create_cif_file_bead_list(luout,tag) 1246 implicit none 1247 integer luout 1248 character*(*) tag 1249 1250#include "rtdb.fh" 1251#include "mafdecls.fh" 1252#include "geom.fh" 1253#include "global.fh" 1254 1255* ***** local variables **** 1256 logical value,isperiodic 1257 integer taglen,i,ii,geom,geomlen,nbeads,nion,i1,i2 1258 real*8 energy,q,rxyz(3),a(3,3),b(3,3),volume 1259 real*8 frac(3),aa,bb,cc,pi,d2,alpha,beta,gmma 1260 1261 character*255 rtdb_name,geom_name 1262 character*2 symbol 1263 character*16 t,name 1264 character*26 dd 1265 1266 integer bead_rtdb 1267 common /bead_list/ bead_rtdb 1268 1269* **** external functions **** 1270 integer inp_strlen,size_bead_list 1271 character*7 bead_index_name 1272 external inp_strlen,size_bead_list 1273 external bead_index_name 1274 1275 1276 taglen = inp_strlen(tag) 1277 nbeads = size_bead_list(tag) 1278 value = geom_create(geom,'bead_tmp_cif') 1279 do i=1,nbeads 1280 1281 rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name' 1282 if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)) 1283 > geom_name = tag(1:taglen)//bead_index_name(i)//':geom' 1284 geomlen = inp_strlen(geom_name) 1285 value = value.and.geom_rtdb_load(bead_rtdb, 1286 > geom,geom_name(1:geomlen)) 1287 value = value.and.geom_ncent(geom,nion) 1288 1289 call amatrix_get_bead_list(tag,i,a) 1290 1291 b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3) 1292 b(2,1) = a(3,2)*a(1,3) - a(1,2)*a(3,3) 1293 b(3,1) = a(1,2)*a(2,3) - a(2,2)*a(1,3) 1294 b(1,2) = a(2,3)*a(3,1) - a(3,3)*a(2,1) 1295 b(2,2) = a(3,3)*a(1,1) - a(1,3)*a(3,1) 1296 b(3,2) = a(1,3)*a(2,1) - a(2,3)*a(1,1) 1297 b(1,3) = a(2,1)*a(3,2) - a(3,1)*a(2,2) 1298 b(2,3) = a(3,1)*a(1,2) - a(1,1)*a(3,2) 1299 b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2) 1300 volume = a(1,1)*b(1,1) 1301 > + a(2,1)*b(2,1) 1302 > + a(3,1)*b(3,1) 1303 1304 volume = 1.0d0/volume 1305 call dscal(9,volume,b,1) 1306 1307* **** determine a,b,c,alpha,beta,gmma *** 1308 pi = 4.0d0*datan(1.0d0) 1309 aa = dsqrt(a(1,1)**2 + a(2,1)**2 +a(3,1)**2) 1310 bb = dsqrt(a(1,2)**2 + a(2,2)**2 +a(3,2)**2) 1311 cc = dsqrt(a(1,3)**2 + a(2,3)**2 +a(3,3)**2) 1312 1313 d2 = (a(1,2)-a(1,3))**2+(a(2,2)-a(2,3))**2+(a(3,2)-a(3,3))**2 1314 alpha = (bb*bb + cc*cc - d2)/(2.0d0*bb*cc) 1315 alpha = dacos(alpha)*180.0d0/pi 1316 1317 d2 = (a(1,3)-a(1,1))**2+(a(2,3)-a(2,1))**2+(a(3,3)-a(3,1))**2 1318 beta = (cc*cc + aa*aa - d2)/(2.0d0*cc*aa) 1319 beta = dacos(beta)*180.0d0/pi 1320 1321 d2 = (a(1,1)-a(1,2))**2+(a(2,1)-a(2,2))**2+(a(3,1)-a(3,2))**2 1322 gmma = (aa*aa + bb*bb - d2)/(2.0d0*aa*bb) 1323 gmma = dacos(gmma)*180.0d0/pi 1324 1325 if (ga_nodeid().eq.0) then 1326 1327 call util_date(dd) 1328 write(luout,1200) 1329 write(luout,1210) dd(1:24) 1330 write(luout,1211) 1331 1332 write(luout,1220) aa * 0.529177d0 1333 write(luout,1221) bb * 0.529177d0 1334 write(luout,1222) cc * 0.529177d0 1335 write(luout,1223) alpha 1336 write(luout,1224) beta 1337 write(luout,1225) gmma 1338 1339 write(luout,1230) 1340 1341 write(luout,1240) 1342 write(luout,1241) 1343 write(luout,1243) 1344 write(luout,1244) 1345 write(luout,1245) 1346 1347 1348 end if 1349 do ii=1,nion 1350 value = value.and.geom_cent_get(geom,ii,t,rxyz,q) 1351 value = value.and. geom_tag_to_element(t,symbol,name,q) 1352 frac(1) = b(1,1)*rxyz(1) 1353 > + b(2,1)*rxyz(2) 1354 > + b(3,1)*rxyz(3) 1355 frac(2) = b(1,2)*rxyz(1) 1356 > + b(2,2)*rxyz(2) 1357 > + b(3,2)*rxyz(3) 1358 frac(3) = b(1,3)*rxyz(1) 1359 > + b(2,3)*rxyz(2) 1360 > + b(3,3)*rxyz(3) 1361 if (ga_nodeid().eq.0) write(luout,1250) symbol,frac 1362 end do 1363 1364 end do 1365 value = value.and.geom_destroy(geom) 1366 1367 1368 return 1369 1200 FORMAT('data_nwchem_pspw') 1370 1210 FORMAT(/'_audit_creation_date ',A) 1371 1211 FORMAT( 1372 > '_audit_creation_method generated by PSPW module of NWChem') 1373 1374 1220 FORMAT(//'_cell_length_a ', F16.4) 1375 1221 FORMAT( '_cell_length_b ', F16.4) 1376 1222 FORMAT( '_cell_length_c ', F16.4) 1377 1223 FORMAT( '_cell_angle_alpha', F16.4) 1378 1224 FORMAT( '_cell_angle_beta ', F16.4) 1379 1225 FORMAT( '_cell_angle_gamma', F16.4) 1380 1381 1230 FORMAT(/'_symmetry_space_group_name_H-M P1 ') 1382 1383 1240 FORMAT(/'loop_') 1384 1241 FORMAT('_atom_site_type_symbol') 1385 1242 FORMAT('_atom_site_label') 1386 1243 FORMAT('_atom_site_fract_x') 1387 1244 FORMAT('_atom_site_fract_y') 1388 1245 FORMAT('_atom_site_fract_z') 1389 1390c 1250 FORMAT(A2,6X,I4,3x,3F14.6) 1391 1250 FORMAT(A2,6X,3F14.6) 1392 end 1393 1394 1395c *********************************************** 1396c * * 1397c * bead_index_name * 1398c * * 1399c *********************************************** 1400 character*7 function bead_index_name(i) 1401 integer i 1402 1403 integer itmp,j0,j1,j2,j3,j4,j5 1404 character*7 name 1405 1406 itmp = i 1407 1408 j5 = itmp/100000 1409 itmp = itmp - j5*100000 1410 j4 = itmp/10000 1411 itmp = itmp - j4*10000 1412 j3 = itmp/1000 1413 itmp = itmp - j3*1000 1414 j2 = itmp/100 1415 itmp = itmp - j2*100 1416 j1 = itmp/10 1417 itmp = itmp - j1*10 1418 j0 = itmp/1 1419 1420 name(1:1) = '_' 1421 name(2:2) = CHAR(j5+ICHAR('0')) 1422 name(3:3) = CHAR(j4+ICHAR('0')) 1423 name(4:4) = CHAR(j3+ICHAR('0')) 1424 name(5:5) = CHAR(j2+ICHAR('0')) 1425 name(6:6) = CHAR(j1+ICHAR('0')) 1426 name(7:7) = CHAR(j0+ICHAR('0')) 1427 bead_index_name = name 1428 return 1429 end 1430 1431c *************************************************** 1432c * * 1433c * bead_includestress * 1434c * * 1435c *************************************************** 1436 logical function bead_includestress() 1437 implicit none 1438 1439#include "rtdb.fh" 1440#include "mafdecls.fh" 1441#include "geom.fh" 1442 1443* ***** local variables **** 1444 logical ostress 1445 1446 integer bead_rtdb 1447 common /bead_list/ bead_rtdb 1448 1449 if (.not.rtdb_get(bead_rtdb,'includestress',mt_log,1,ostress)) 1450 > ostress = .false. 1451 1452 bead_includestress = ostress 1453 return 1454 end 1455 1456