1c 2c $Id$ 3c 4 5 6* ***************************************************** 7* * * 8* * pspw_dplot_input * 9* * * 10* ***************************************************** 11 12 subroutine pspw_dplot_input(rtdb) 13 implicit none 14 integer rtdb 15c 16#include "errquit.fh" 17#include "inp.fh" 18#include "bafdecls.fh" 19#include "rtdb.fh" 20#include "nwc_const.fh" 21c 22 integer num_dirs ! No. of known directives 23 parameter (num_dirs = 16) 24 character*22 dirs(num_dirs) 25 data dirs / 'vectors', 26 > 'density', 27 > 'orbital', 28 > 'position_tolerance', 29 > 'elf', 30 > '2d_grid', 31 > '3d_grid', 32 > 'origin', 33 > 'limitxyz', 34 > 'ncell', 35 > 'atom_truncate', 36 > '1d_grid', 37 > 'orbital2', 38 > 'density_matrix', 39 > 'structure_factor', 40 > 'end'/ 41 42 integer num_dnames ! No. of density directives 43 parameter (num_dnames = 8) 44 character*22 dnames(num_dnames) 45 data dnames / 'total', 46 > 'diff', 47 > 'alpha', 48 > 'beta', 49 > 'laplacian', 50 > 'potential', 51 > 'up', 52 > 'down' / 53 integer num_enames ! No. of ELF directives 54 parameter (num_enames = 5) 55 character*22 enames(num_enames) 56 data enames / 'restricted', 57 > 'up', 58 > 'down', 59 > 'alpha', 60 > 'beta' / 61 62 logical value 63 integer ind ! Index of matched directive 64 integer number,number2,count,ia 65 integer name1_len,name2_len,name3_len,name4_len 66 character*50 name1,name2,name3,name4 67 character*50 filename 68 character*50 wavefunction_filename 69 character*255 test 70 real*8 position_tolerance 71 real*8 o(3),x(3),y(3),z(3) 72 real*8 sizex(2),sizey(2),sizez(2),scal 73 integer j,jstart,jlast,jstride 74 integer nx,ny,nz,ncell(3),idx(2) 75c 76c 77* *** initializations **** 78 position_tolerance=0.001d0 79 call util_file_prefix('movecs',wavefunction_filename) 80 ia = ICHAR('a') 81 count = 0 82 value = rtdb_put(rtdb,'pspw_dplot:count',mt_int,1,count) 83 if (.not.value) 84 > call errquit( 85 > 'pspw_dplot_input: rtdb_put failed for count', 0, RTDB_ERR) 86 87 88 10 if (.not. inp_read()) 89 > call errquit( 90 > 'pspw_dplot_input: inp_read failed', 0, INPUT_ERR) 91 if (.not. inp_a(test)) 92 > call errquit( 93 > 'pspw_dplot_input: failed to read keyword', 0, 94 & INPUT_ERR) 95 if (.not. inp_match(num_dirs, .false., test, dirs, ind)) 96 > call errquit( 97 > 'pspw_dplot_input: unknown directive', 0, INPUT_ERR) 98 99 100 goto ( 100, 200, 300, 400, 500, 600,700,800,900,1000,1100,1200, 101 > 350, 360, 370, 102 > 9999) ind 103 call errquit( 104 > 'pspw_dplot_input: unimplemented directive', ind, INPUT_ERR) 105 106 107c 108c vectors 109c 110 100 if (.not. inp_a(wavefunction_filename)) 111 > call errquit( 112 > 'pspw_dplot_input: failed to read vector', 0, INPUT_ERR) 113 goto 10 114 115c 116c density 117c 118* *** read density option **** 119 200 if (.not. inp_a(test)) 120 > call errquit( 121 > 'pspw_dplot_input: failed to read density keyword', 0, 122 & INPUT_ERR) 123 124* *** density number **** 125 if (.not.inp_match(num_dnames,.false.,test,dnames,number)) 126 > number = 1 127 if (number .eq. 7) number = 3 128 if (number .eq. 8) number = 4 129* !*** number = 1 - total 130* !*** number = 2 - difference 131* !*** number = 3 - alpha 132* !*** number = 4 - beta 133* !*** number = 5 - laplacian 134* !*** number = 6 - potential 135* !*** number = 7 - ELF 136 137* **** make density numbers negative **** 138 number = -number 139 140* *** read filename **** 141 if (.not. inp_a(filename)) 142 > call errquit( 143 > 'pspw_dplot_input: failed to read density filename', 0, 144 & INPUT_ERR) 145 146* **** define name - not very elegent and could break if **** 147* **** count becomes very large **** 148 count = count + 1 149 name1 = 'pspw_dplot:filename'//CHAR(count-1+ia) 150 name2 = 'pspw_dplot:number'//CHAR(count-1+ia) 151 name1_len = index(name1,' ') - 1 152 name2_len = index(name2,' ') - 1 153 154 ind = index(filename,' ') - 1 155 value = rtdb_cput(rtdb,name1(1:name1_len),1,filename(1:ind)) 156 value = value.and.rtdb_put(rtdb,name2(1:name2_len), 157 > mt_int,1,number) 158 value = value.and.rtdb_put(rtdb,'pspw_dplot:count', 159 > mt_int,1,count) 160 if (.not.value) 161 > call errquit( 162 > 'pspw_dplot_input: rtdb_put failed for density', 0, 163 & RTDB_ERR) 164 165 goto 10 166 167 168c 169c orbital 170c 171* *** read orbital number **** 172 300 if (.not. inp_i(number)) 173 > call errquit( 174 > 'pspw_dplot_input: failed to read orbital number', 0, 175 & INPUT_ERR) 176 177* *** read filename **** 178 if (.not. inp_a(filename)) 179 > call errquit( 180 > 'pspw_dplot_input: failed to read orbital filename', 0, 181 & INPUT_ERR) 182 183* **** define name - not very elegent and could break if **** 184* **** count becomes very large **** 185 count = count + 1 186 name1 = 'pspw_dplot:filename'//CHAR(count-1+ia) 187 name2 = 'pspw_dplot:number'//CHAR(count-1+ia) 188 name1_len = index(name1,' ') - 1 189 name2_len = index(name2,' ') - 1 190 191 ind = index(filename,' ') - 1 192 value = rtdb_cput(rtdb,name1(1:name1_len),1,filename(1:ind)) 193 value = value.and.rtdb_put(rtdb,name2(1:name2_len), 194 > mt_int,1,number) 195 value = value.and.rtdb_put(rtdb,'pspw_dplot:count', 196 > mt_int,1,count) 197 if (.not.value) 198 > call errquit( 199 > 'pspw_dplot_input: rtdb_put failed for orbital', 0, RTDB_ERR) 200 201 goto 10 202 203 204c 205c orbital2 206c 207* *** read orbital number **** 208 350 if (.not. inp_i(number)) 209 > call errquit( 210 > 'pspw_dplot_input: failed to read orbital number1',0, 211 & INPUT_ERR) 212 if (.not. inp_i(number2)) 213 > call errquit( 214 > 'pspw_dplot_input: failed to read orbital number2',0, 215 & INPUT_ERR) 216 217* *** read filename **** 218 if (.not. inp_a(filename)) 219 > call errquit( 220 > 'pspw_dplot_input: failed to read orbital2 filename',0, 221 & INPUT_ERR) 222 223* **** define name - not very elegent and could break if **** 224* **** count becomes very large **** 225 count = count + 1 226 name1 = 'pspw_dplot:filename'//CHAR(count-1+ia) 227 name2 = 'pspw_dplot:number1'//CHAR(count-1+ia) 228 name3 = 'pspw_dplot:number2'//CHAR(count-1+ia) 229 name1_len = index(name1,' ') - 1 230 name2_len = index(name2,' ') - 1 231 name3_len = index(name3,' ') - 1 232 233 ind = index(filename,' ') - 1 234 value = rtdb_cput(rtdb,name1(1:name1_len),1,filename(1:ind)) 235 value = value.and.rtdb_put(rtdb,name2(1:name2_len), 236 > mt_int,1,number) 237 value = value.and.rtdb_put(rtdb,name3(1:name3_len), 238 > mt_int,1,number2) 239 value = value.and.rtdb_put(rtdb,'pspw_dplot:count', 240 > mt_int,1,count) 241 if (.not.value) 242 > call errquit( 243 > 'pspw_dplot_input:rtdb_put failed for orbital2',0,RTDB_ERR) 244 245 goto 10 246 247c 248c density_matrix ms x,y,z 249c 250* *** read ms,x,y,z **** 251 360 if (.not. inp_i(number2)) 252 > call errquit( 253 > 'pspw_dplot_input: failed to read density_matrix ms', 254 > 0,INPUT_ERR) 255 if (.not. inp_f(x(1))) 256 > call errquit('pspw_dplot_input:failed to read density_matrix x', 257 > 0,INPUT_ERR) 258 if (.not. inp_f(x(2))) 259 > call errquit('pspw_dplot_input:failed to read density_matrix y', 260 > 0,INPUT_ERR) 261 if (.not. inp_f(x(3))) 262 > call errquit('pspw_dplot_input:failed to read density_matrix z', 263 > 0,INPUT_ERR) 264 265* *** read filename **** 266 if (.not. inp_a(filename)) 267 > call errquit( 268 > 'pspw_dplot_input: failed to read density_matrix filename', 269 > 0,INPUT_ERR) 270 271* **** define name - not very elegent and could break if **** 272* **** count becomes very large **** 273 count = count + 1 274 name1 = 'pspw_dplot:filename'//CHAR(count-1+ia) 275 name2 = 'pspw_dplot:number'//CHAR(count-1+ia) 276 name3 = 'pspw_dplot:dmatrix_ms'//CHAR(count-1+ia) 277 name4 = 'pspw_dplot:dmatrix_xyz'//CHAR(count-1+ia) 278 name1_len = index(name1,' ') - 1 279 name2_len = index(name2,' ') - 1 280 name3_len = index(name3,' ') - 1 281 name4_len = index(name4,' ') - 1 282 283 number = -9 284 ind = index(filename,' ') - 1 285 value = rtdb_cput(rtdb,name1(1:name1_len),1,filename(1:ind)) 286 value = value.and.rtdb_put(rtdb,name2(1:name2_len), 287 > mt_int,1,number) 288 value = value.and.rtdb_put(rtdb,name3(1:name3_len), 289 > mt_int,1,number2) 290 value = value.and.rtdb_put(rtdb,name4(1:name4_len), 291 > mt_dbl,3,x) 292 value = value.and.rtdb_put(rtdb,'pspw_dplot:count', 293 > mt_int,1,count) 294 if (.not.value) 295 > call errquit( 296 > 'pspw_dplot_input:rtdb_put failed for densty_dmatix',0,RTDB_ERR) 297 298 goto 10 299 300 301 302c 303c structure_factor 304c 305* *** read structure factor option **** 306 370 if (.not. inp_a(filename)) 307 > call errquit( 308 > 'pspw_dplot_input: failed read structure_factor filename', 0, 309 & INPUT_ERR) 310 311 number = -10 312* **** define name - not very elegent and could break if **** 313* **** count becomes very large **** 314 count = count + 1 315 name1 = 'pspw_dplot:filename'//CHAR(count-1+ia) 316 name2 = 'pspw_dplot:number'//CHAR(count-1+ia) 317 name1_len = index(name1,' ') - 1 318 name2_len = index(name2,' ') - 1 319 320 ind = index(filename,' ') - 1 321 value = rtdb_cput(rtdb,name1(1:name1_len),1,filename(1:ind)) 322 value = value.and.rtdb_put(rtdb,name2(1:name2_len), 323 > mt_int,1,number) 324 value = value.and.rtdb_put(rtdb,'pspw_dplot:count', 325 > mt_int,1,count) 326 if (.not.value) 327 > call errquit( 328 > 'pspw_dplot_input: rtdb_put failed for structure_factor', 0, 329 & RTDB_ERR) 330 331 goto 10 332 333 334 335 336c 337c position_tolerance 338c 339* **** read position_tolerance **** 340 400 if (.not. inp_f(position_tolerance)) 341 > call errquit( 342 > 'pspw_dplot_input: failed to read position tolerance', 0, 343 & INPUT_ERR) 344 345 goto 10 346 347c 348c ELF 349c 350 500 if (.not. inp_a(test)) 351 > call errquit( 352 > 'pspw_dplot_input: failed to read ELF keyword', 0, 353 & INPUT_ERR) 354 355* *** ELF number **** 356 if (.not.inp_match(num_enames,.false.,test,enames,number)) 357 > number = -7 358 if (number .eq. 1) number = -7 359 if (number .eq. 2) number = -7 360 if (number .eq. 3) number = -8 361 if (number .eq. 4) number = -7 362 if (number .eq. 5) number = -8 363* !*** number = -7 - ELF up 364* !*** number = -8 - ELF down 365 366* *** read ELF filename **** 367 if (.not. inp_a(filename)) 368 > filename = test 369 370* **** define name - not very elegent and could break if **** 371* **** count becomes very large **** 372 count = count + 1 373 name1 = 'pspw_dplot:filename'//CHAR(count-1+ia) 374 name2 = 'pspw_dplot:number'//CHAR(count-1+ia) 375 name1_len = index(name1,' ') - 1 376 name2_len = index(name2,' ') - 1 377 378 ind = index(filename,' ') - 1 379 value = rtdb_cput(rtdb,name1(1:name1_len),1,filename(1:ind)) 380 value = value.and.rtdb_put(rtdb,name2(1:name2_len), 381 > mt_int,1,number) 382 value = value.and.rtdb_put(rtdb,'pspw_dplot:count', 383 > mt_int,1,count) 384 if (.not.value) 385 > call errquit( 386 > 'pspw_dplot_input: rtdb_put failed for ELF', 0, RTDB_ERR) 387 388 389 goto 10 390 391c 392c 2d_grid 393c 394* **** read o,x,sizex,y,sizey **** 395 600 value = inp_read() 396 value = value.and.inp_f(o(1)) 397 value = value.and.inp_f(o(2)) 398 value = value.and.inp_f(o(3)) 399 400 value = value.and.inp_read() 401 value = value.and.inp_f(x(1)) 402 value = value.and.inp_f(x(2)) 403 value = value.and.inp_f(x(3)) 404 value = value.and.inp_f(sizex(1)) 405 value = value.and.inp_f(sizex(2)) 406 407 value = value.and.inp_read() 408 value = value.and.inp_f(y(1)) 409 value = value.and.inp_f(y(2)) 410 value = value.and.inp_f(y(3)) 411 value = value.and.inp_f(sizey(1)) 412 value = value.and.inp_f(sizey(2)) 413 414 415 value = value.and.rtdb_put(rtdb,'pspw_dplot:2d_grid:o', 416 > mt_dbl,3,o) 417 value = value.and.rtdb_put(rtdb,'pspw_dplot:2d_grid:x', 418 > mt_dbl,3,x) 419 value = value.and.rtdb_put(rtdb,'pspw_dplot:2d_grid:y', 420 > mt_dbl,3,y) 421 value = value.and. 422 > rtdb_put(rtdb,'pspw_dplot:2d_grid:sizex', 423 > mt_dbl,2,sizex) 424 value = value.and. 425 > rtdb_put(rtdb,'pspw_dplot:2d_grid:sizey', 426 > mt_dbl,2,sizey) 427 428 if (.not. value) 429 > call errquit( 430 > 'pspw_dplot_input: 2d_grid failed to read', 0, RTDB_ERR) 431 432 goto 10 433 434c 435c 3d_grid 436c 437* **** read o,x,sizex,y,sizey,y,sizez **** 438 700 value = inp_read() 439 value = value.and.inp_f(o(1)) 440 value = value.and.inp_f(o(2)) 441 value = value.and.inp_f(o(3)) 442 443 value = value.and.inp_read() 444 value = value.and.inp_f(x(1)) 445 value = value.and.inp_f(x(2)) 446 value = value.and.inp_f(x(3)) 447 value = value.and.inp_f(sizex(1)) 448 value = value.and.inp_f(sizex(2)) 449 value = value.and.inp_i(nx) 450 451 value = value.and.inp_read() 452 value = value.and.inp_f(y(1)) 453 value = value.and.inp_f(y(2)) 454 value = value.and.inp_f(y(3)) 455 value = value.and.inp_f(sizey(1)) 456 value = value.and.inp_f(sizey(2)) 457 value = value.and.inp_i(ny) 458 459 value = value.and.inp_read() 460 value = value.and.inp_f(z(1)) 461 value = value.and.inp_f(z(2)) 462 value = value.and.inp_f(z(3)) 463 value = value.and.inp_f(sizez(1)) 464 value = value.and.inp_f(sizez(2)) 465 value = value.and.inp_i(nz) 466 467 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:o', 468 > mt_dbl,3,o) 469 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:x', 470 > mt_dbl,3,x) 471 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:y', 472 > mt_dbl,3,y) 473 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:z', 474 > mt_dbl,3,z) 475 value = value.and. 476 > rtdb_put(rtdb,'pspw_dplot:3d_grid:sizex', 477 > mt_dbl,2,sizex) 478 value = value.and. 479 > rtdb_put(rtdb,'pspw_dplot:3d_grid:sizey', 480 > mt_dbl,2,sizey) 481 value = value.and. 482 > rtdb_put(rtdb,'pspw_dplot:3d_grid:sizez', 483 > mt_dbl,2,sizez) 484 value = value.and. 485 > rtdb_put(rtdb,'pspw_dplot:3d_grid:nx', 486 > mt_int,1,nx) 487 value = value.and. 488 > rtdb_put(rtdb,'pspw_dplot:3d_grid:ny', 489 > mt_int,1,ny) 490 value = value.and. 491 > rtdb_put(rtdb,'pspw_dplot:3d_grid:nz', 492 > mt_int,1,nz) 493 494 if (.not. value) 495 > call errquit( 496 > 'pspw_dplot_input: 3d_grid failed to read', 0, RTDB_ERR) 497 498 goto 10 499 500 501c 502c translate_origin 503c 504 800 value = inp_f(o(1)) 505 value = value.and.inp_f(o(2)) 506 value = value.and.inp_f(o(3)) 507 508 value = value.and.rtdb_put(rtdb,'pspw_dplot:origin', 509 > mt_dbl,3,o) 510 511 if (.not. value) 512 > call errquit( 513 > 'pspw_dplot_input: translate_origin failed to read',0,0) 514 515 goto 10 516 517 518c 519c limitxyz 520c 521* **** read o,x,sizex,y,sizey,y,sizez **** 522 900 call get_scalefrominput(scal) 523 524 value = value.and.inp_read() 525 value = value.and.inp_f(sizex(1)) 526 value = value.and.inp_f(sizex(2)) 527 value = value.and.inp_i(nx) 528 sizex(1) = scal*sizex(1) 529 sizex(2) = scal*sizex(2) 530 531 value = value.and.inp_read() 532 value = value.and.inp_f(sizey(1)) 533 value = value.and.inp_f(sizey(2)) 534 value = value.and.inp_i(ny) 535 sizey(1) = scal*sizey(1) 536 sizey(2) = scal*sizey(2) 537 538 value = value.and.inp_read() 539 value = value.and.inp_f(sizez(1)) 540 value = value.and.inp_f(sizez(2)) 541 value = value.and.inp_i(nz) 542 sizez(1) = scal*sizez(1) 543 sizez(2) = scal*sizez(2) 544 545 !*** set origin and axes *** 546 o(1) = 0.0d0 547 o(2) = 0.0d0 548 o(3) = 0.0d0 549 550 x(1) = 1.0d0 551 x(2) = 0.0d0 552 x(3) = 0.0d0 553 554 y(1) = 0.0d0 555 y(2) = 1.0d0 556 y(3) = 0.0d0 557 558 z(1) = 0.0d0 559 z(2) = 0.0d0 560 z(3) = 1.0d0 561 562 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:o', 563 > mt_dbl,3,o) 564 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:x', 565 > mt_dbl,3,x) 566 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:y', 567 > mt_dbl,3,y) 568 value = value.and.rtdb_put(rtdb,'pspw_dplot:3d_grid:z', 569 > mt_dbl,3,z) 570 value = value.and. 571 > rtdb_put(rtdb,'pspw_dplot:3d_grid:sizex', 572 > mt_dbl,2,sizex) 573 value = value.and. 574 > rtdb_put(rtdb,'pspw_dplot:3d_grid:sizey', 575 > mt_dbl,2,sizey) 576 value = value.and. 577 > rtdb_put(rtdb,'pspw_dplot:3d_grid:sizez', 578 > mt_dbl,2,sizez) 579 value = value.and. 580 > rtdb_put(rtdb,'pspw_dplot:3d_grid:nx', 581 > mt_int,1,nx) 582 value = value.and. 583 > rtdb_put(rtdb,'pspw_dplot:3d_grid:ny', 584 > mt_int,1,ny) 585 value = value.and. 586 > rtdb_put(rtdb,'pspw_dplot:3d_grid:nz', 587 > mt_int,1,nz) 588 589 if (.not. value) 590 > call errquit( 591 > 'pspw_dplot_input: failed to write limitxyz',0,RTDB_ERR) 592 593 goto 10 594 595 596c 597c ncell 598c 599 1000 if (.not.inp_i(ncell(1))) ncell(1) = 0 600 if (.not.inp_i(ncell(2))) ncell(2) = 0 601 if (.not.inp_i(ncell(3))) ncell(3) = 0 602 if (.not.rtdb_put(rtdb,'pspw_dplot:ncell',mt_int,3,ncell)) 603 > call errquit( 604 > 'pspw_dplot_input: failed to write ncell',0,RTDB_ERR) 605 606 goto 10 607 608c 609c atom_truncate 610c 611 1100 if (.not.BA_push_get(mt_int,nw_max_atom,'idx',idx(2),idx(1))) 612 > call errquit( 613 > 'pspw_dplot_input:failed allocating idx',0,MA_ERR) 614 615 nx = 0 616 do while (inp_irange(jstart,jlast,jstride)) 617 do j=jstart,jlast,jstride 618 int_mb(idx(1)+nx) = j 619 nx = nx+1 620 end do 621 end do 622 if (.not.rtdb_put(rtdb,'pspw_dplot:atom_truncate_size', 623 > mt_int,1,nx)) 624 > call errquit('pspw_dplot_input:failed write atom_truncate_size', 625 > 0,RTDB_ERR) 626 if (.not.rtdb_put(rtdb,'pspw_dplot:atom_truncate', 627 > mt_int,nx,int_mb(idx(1)))) 628 > call errquit('pspw_dplot_input:failed writing atom_truncate', 629 > 0,RTDB_ERR) 630 if (.not.BA_pop_stack(idx(2))) 631 > call errquit('pspw_dplot_input:failed deallocating idx', 632 > 0,MA_ERR) 633 634 goto 10 635 636c 637c 1d_grid 638c 639* **** read o,x,sizex,y,sizey **** 640 1200 value = inp_read() 641 value = value.and.inp_i(nx) 642 643 value = value.and.inp_read() 644 value = value.and.inp_f(o(1)) 645 value = value.and.inp_f(o(2)) 646 value = value.and.inp_f(o(3)) 647 648 value = value.and.inp_read() 649 value = value.and.inp_f(x(1)) 650 value = value.and.inp_f(x(2)) 651 value = value.and.inp_f(x(3)) 652 653 654 value = value.and.rtdb_put(rtdb,'pspw_dplot:1d_grid:o', 655 > mt_dbl,3,o) 656 value = value.and.rtdb_put(rtdb,'pspw_dplot:1d_grid:x', 657 > mt_dbl,3,x) 658 value = value.and. 659 > rtdb_put(rtdb,'pspw_dplot:1d_grid:nx', 660 > mt_int,1,nx) 661 662 663 goto 10 664 665 666 9999 continue 667 668 ind = index(wavefunction_filename,' ') - 1 669 value = rtdb_cput(rtdb,'pspw_dplot:wavefunction_filename', 670 > 1,wavefunction_filename(1:ind)) 671 value = value.and. 672 > rtdb_put(rtdb,'pspw_dplot:position_tolerance', 673 > mt_dbl,1,position_tolerance) 674 if (.not.value) 675 > call errquit( 676 > 'pspw_dplot_input: rtdb_put failed for vector', 0, RTDB_ERR) 677 678 return 679 end 680