1C 2C rt_tddft_input_field.F 3C 4C Parses input deck for rt-tddft field (excitation) parameters. 5C 6C 7 subroutine rt_tddft_input_field (rtdb, field_name, nfields) 8 implicit none 9 10#include "rt_tddft.fh" 11#include "errquit.fh" 12#include "inp.fh" 13#include "rtdb.fh" 14#include "mafdecls.fh" 15#include "stdio.fh" 16 17 18C == Inputs == 19 integer, intent(in) :: rtdb 20 character*16, intent(in) :: field_name !hardcoded to match geom name max size 21 integer, intent(in) :: nfields !this is the number of the current field 22 23 24C == Parameters == 25 character(*), parameter :: pname = "rt_tddft_input_field: " 26 27 28C == Variables == 29 logical done 30 character*255 test, filename 31 32 type (rt_field_t) prev_field, this_field 33 integer i 34 35 character*20 type 36 character*20 spin 37 character spin1 38 double precision max 39 double precision center 40 double precision frequency 41 double precision width 42 double precision phase 43 double precision start 44 character*16 polarization !x,y,z for dipole; xx,xy,xz,... for quad 45 double precision theta, phi 46 47 logical lhave_center 48 logical lhave_start 49 logical lhave_polarization 50 logical lhave_width 51 logical lhave_max 52 logical lhave_type 53 logical lhave_frequency 54 logical lhave_spin 55 logical lhave_phase 56 57 if (nfields .gt. rt_max_fields) 58 $ call errquit (pname//"cannot exceed max num fields", 0, 0) 59 60 61 lhave_center = .false. 62 lhave_polarization = .false. 63 lhave_width = .false. 64 lhave_type = .false. 65 lhave_frequency = .false. 66 lhave_spin = .false. 67 lhave_phase = .false. 68 lhave_max = .false. 69 lhave_start = .false. 70 71 phase = 0d0 ! default phase = 0 72 73C 74C Parse the input; we will put in rtdb later after checking. 75C 76 done = .false. 77 do while (.not. done) 78 79 if (.not. inp_read()) 80 $ call errquit(pname//'Read failed input',0, INPUT_ERR) 81 if (.not. inp_a(test)) 82 $ call errquit(pname//'Read failed keyword',0, INPUT_ERR) 83 84 85 86C 87C type (delta, cw, gaussian) 88C 89 if (inp_compare(.false.,test,'type')) then 90 if (.not. inp_a (type)) 91 $ call errquit (pname//"failed to read field type",0,0) 92 93 if ( type .eq. "file") then 94 if (.not. inp_a (filename)) call errquit (pname// 95 $ "failed to read field file name", 0, 0) 96 97 else 98 if ( (type.ne."cw").and. 99 $ (type.ne."delta").and. 100 $ (type.ne."hann").and. 101 $ (type.ne."sin2ramp").and. 102 $ (type.ne."gaussian")) 103 $ call errquit (pname//"invalid field type: "// 104 $ type,0,0) 105 endif 106 lhave_type = .true. 107 108 109C 110C spin which the field acts on 111C 112 elseif (inp_compare(.false.,test,'spin')) then 113 if (.not. inp_a (spin)) 114 $ call errquit (pname// 115 $ "failed to read field target spin",0,0) 116 117 lhave_spin = .true. 118 119 120 121C 122C max value of the field 123C 124 elseif (inp_compare(.false.,test,'max')) then 125 if (.not.inp_f(max)) call errquit (pname// 126 $ "max takes a float", 0, 0) 127 lhave_max = .true. 128 129C 130C center the field (only for gaussian and Hann) 131C 132 elseif (inp_compare(.false.,test,'center')) then 133 if (.not.inp_f(center)) call errquit (pname// 134 $ "center takes a float >= 0", 0, 0) 135 lhave_center = .true. 136 137 138C 139C start of the ramp up (only for ramps) 140C 141 elseif (inp_compare(.false.,test,'start')) then 142 if (.not.inp_f(start)) call errquit (pname// 143 $ "start takes a float >= 0", 0, 0) 144 lhave_start = .true. 145 146 147C 148C width the field (only for gaussian and Hann) 149C 150 elseif (inp_compare(.false.,test,'width')) then 151 if (.not.inp_f(width)) call errquit (pname// 152 $ "width takes a float >= 0", 0, 0) 153 lhave_width = .true. 154 155 156C 157C frequency the field (only for gaussian and cw) 158C 159 elseif (inp_compare(.false.,test,'frequency')) then 160 if (.not.inp_f(frequency)) call errquit (pname// 161 $ "frequency takes a float >= 0", 0, 0) 162 lhave_frequency = .true. 163 164 165C 166C field polarization 167C 168c$$$ elseif (inp_compare(.false.,test,'polarization')) then 169c$$$ if (.not.inp_a(polarization)) call errquit (pname// 170c$$$ $ "polarization can be: x,y,z (for dipole); "// 171c$$$ $ "xx,xy,xz,... (for quad)", 0, 0) 172c$$$ lhave_polarization = .true. 173 174 elseif (inp_compare(.false.,test,'polarization')) then 175 if (.not.inp_a(polarization)) call errquit (pname// 176 $ "polarization can be: x,y,z, angle", 0, 0) 177 178 if (inp_compare(.false.,polarization,'angle')) then !need to get theta, phi values 179 if (.not.inp_f(theta)) call errquit (pname// 180 $ "angle theta takes a float", 0, 0) 181 if (.not.inp_f(phi)) call errquit (pname// 182 $ "angle phi takes a float", 0, 0) 183 endif 184 lhave_polarization = .true. 185 186 187C 188C phase (only for gaussian and cw) 189C 190 elseif (inp_compare(.false.,test,'phase')) then 191 if (.not.inp_f(phase)) call errquit (pname// 192 $ "phase takes a float >= 0", 0, 0) 193 lhave_phase = .true. 194 195 196C 197C end of parse 198C 199 else if (inp_compare(.false.,test,'end')) then 200 done = .true. 201 else 202 call errquit(pname//'Unknown directive: '//trim(test), 203 $ 0, INPUT_ERR) 204 endif 205 206 enddo 207 208 209C 210C Now check that we have all required parameters, no superfluous 211C ones, no name clashes with other fields, and that params are 212C reasonable (e.g., no negative times, etc). 213C 214 if (nfields .gt. 1) then 215 do i = 1, nfields - 1 216 call rt_tddft_field_rtdb_get (rtdb, i, prev_field) 217 if (prev_field%name .eq. field_name) 218 $ call errquit (pname//"cannot have multiple fields"// 219 $ " with the same name: "//trim(field_name), 0, 0) 220 enddo 221 endif 222 223 if (.not. lhave_type) 224 $ call errquit (pname//trim(field_name)// 225 $ ": must supply a field type", 0, 0) 226 227 if (lhave_spin) then 228 if (spin.eq."alpha") then 229 spin1 = "a" 230 elseif (spin.eq."beta") then 231 spin1 = "b" 232 elseif (spin.eq."total") then 233 spin1 = "t" 234 else 235 spin1 = "X" 236 call errquit (pname//"invalid field spin: "//spin,0,0) 237 endif 238 else 239 spin1 = "t" !default to acting on all spins endif 240 endif 241 242 if (.not. lhave_spin) spin = "total" !default to acting on both spins 243 244 245 if ( type .eq. "file" ) then 246 if (lhave_max) call errquit (pname//trim(field_name)// 247 $ ": cannot specify field max if reading from file", 0, 0) 248 if (lhave_polarization) call errquit (pname//trim(field_name)// 249 $ ": cannot specify polarization if reading from file", 0,0) 250 if (lhave_frequency) call errquit (pname//trim(field_name)// 251 $ ": cannot specify frequency if reading from file", 0, 0) 252 if (lhave_center) call errquit (pname//trim(field_name)// 253 $ ": cannot specify center if reading from file", 0, 0) 254 if (lhave_width) call errquit (pname//trim(field_name)// 255 $ ": cannot specify width if reading from file", 0, 0) 256 if (lhave_phase) call errquit (pname//trim(field_name)// 257 $ ": cannot specify phase if reading from file", 0, 0) 258 259 else ! internal (non-file) fields 260 if (.not. lhave_max) 261 $ call errquit (pname//trim(field_name)// 262 $ ": must supply a field max", 0, 0) 263 264 if (.not. lhave_polarization) 265 $ call errquit (pname//trim(field_name)// 266 $ ": must supply a field polarization", 0,0) 267 268 if (type .eq. "cw") then 269 if (.not. lhave_frequency) 270 $ call errquit (pname//trim(field_name)// 271 $ ": must supply a frequency if doing cw", 0,0) 272 273 if (lhave_center) call errquit (pname//trim(field_name)// 274 $ ": cannot specify center if cw", 0,0) 275 276 if (lhave_width) call errquit (pname//trim(field_name)// 277 $ ": cannot specify width if cw", 0,0) 278 279 if (lhave_start) call errquit (pname// 280 $ trim(field_name)// 281 $ ": cannot specify start if doing cw", 0, 0) 282 endif 283 284 if (type .eq. "gaussian") then 285 if (.not. lhave_frequency) 286 $ call errquit (pname//trim(field_name)// 287 $ ": must supply a frequency if doing gaussian", 0,0) 288 289 if (.not. lhave_center) call errquit (pname// 290 $ trim(field_name)// 291 $ ": must specify center if gaussian", 0,0) 292 293 if (.not. lhave_width) call errquit (pname// 294 $ trim(field_name)// 295 $ ": must specify width if gaussian", 0,0) 296 297 if (lhave_start) call errquit (pname// 298 $ trim(field_name)// 299 $ ": cannot specify start if doing gaussian", 0, 0) 300 endif 301 302 if (type .eq. "hann") then 303 if (.not. lhave_frequency) 304 $ call errquit (pname//trim(field_name)// 305 $ ": must supply a frequency if doing Hann", 0,0) 306 307 if (.not. lhave_center) call errquit (pname// 308 $ trim(field_name)// 309 $ ": must specify center if Hann", 0,0) 310 311 if (.not. lhave_width) call errquit (pname// 312 $ trim(field_name)// 313 $ ": must specify width if Hann", 0,0) 314 315 if (lhave_start) call errquit (pname// 316 $ trim(field_name)// 317 $ ": cannot specify start if doing hann", 0, 0) 318 endif 319 320 321 if (type .eq. "delta") then 322 if (lhave_frequency) call errquit (pname//trim(field_name)// 323 $ ": cannot supply a frequency if doing delta", 0,0) 324 325 if (.not. lhave_center) then 326 center = 0d0 !default delta kick to t=0 327 lhave_center = .true. 328 endif 329 330 if (lhave_width) call errquit (pname//trim(field_name)// 331 $ ": cannot specify width if delta", 0,0) 332 333 if (lhave_start) call errquit (pname// 334 $ trim(field_name)// 335 $ ": cannot specify start if doing delta", 0, 0) 336 337 endif 338 339 340 if (type .eq. "sin2ramp") then 341 if (.not. lhave_frequency) 342 $ call errquit (pname//trim(field_name)// 343 $ ": must supply a frequency if doing sin2ramp", 0,0) 344 345 if (.not. lhave_width) call errquit (pname// 346 $ trim(field_name)// 347 $ ": must specify width if doing sin2ramp", 0,0) 348 349 if (.not. lhave_start) call errquit (pname// 350 $ trim(field_name)// 351 $ ": must specify start if doing sin2ramp", 0, 0) 352 353 if (lhave_center) call errquit (pname// 354 $ trim(field_name)// 355 $ ": cannot specify center if doing sin2ramp", 0, 0) 356 endif 357 358 359 if ( (polarization.ne."x").and. 360 $ (polarization.ne."y").and. 361 $ (polarization.ne."z").and. 362 $ (polarization.ne."angle")) 363 $ call errquit (pname//trim(field_name)// 364 $ ": polarization must be x, y, z, or angle", 365 $ 0,0) 366 367 if ( (lhave_center).and.(center.lt.0d0) ) 368 $ call errquit (pname//trim(field_name)// 369 $ ": center must be positive", 0, 0) 370 371 if ( (lhave_width).and.(width.lt.0d0) ) 372 $ call errquit (pname//trim(field_name)// 373 $ ": width must be positive", 0, 0) 374 375C 376C Frequency-related stuff only valid for CW and pulses (gaussian, hann) 377C 378 if (lhave_phase .or. lhave_frequency) then 379 if ((type .ne. "cw").and.(type .ne. "gaussian") 380 $ .and. (type .ne. "hann").and.(type .ne. "sin2ramp")) 381 $ call errquit (pname// 382 $ "phase and frequency only valid for "// 383 $ "CW, gaussian, hann, and sin2ramp",0,0) 384 endif 385 endif 386 387C 388C Load into rtdb 389C 390 this_field%name = field_name 391 this_field%type = type 392 393 if (type.eq."file") then 394 this_field%filename = filename 395 elseif (type.eq."cw") then 396 this_field%polarization = polarization 397 this_field%max = max 398 this_field%spin = spin1 399 this_field%frequency = frequency 400 this_field%phase = phase 401 this_field%width = -99d0 402 this_field%center = -99d0 403 this_field%start = -99d0 404 elseif (type.eq."gaussian") then 405 this_field%polarization = polarization 406 this_field%max = max 407 this_field%spin = spin1 408 this_field%frequency = frequency 409 this_field%phase = phase 410 this_field%width = width 411 this_field%center = center 412 this_field%start = -99d0 413 elseif (type.eq."hann") then 414 this_field%polarization = polarization 415 this_field%max = max 416 this_field%spin = spin1 417 this_field%frequency = frequency 418 this_field%phase = phase 419 this_field%width = width 420 this_field%center = center 421 this_field%start = -99d0 422 elseif (type.eq."sin2ramp") then 423 this_field%polarization = polarization 424 this_field%max = max 425 this_field%spin = spin1 426 this_field%frequency = frequency 427 this_field%phase = phase 428 this_field%width = width 429 this_field%start = start 430 this_field%center = -99d0 431 elseif (type.eq."delta") then 432 this_field%polarization = polarization 433 this_field%max = max 434 this_field%spin = spin1 435 this_field%frequency = -99d0 436 this_field%width = -99d0 437 this_field%center = center 438 this_field%start = -99d0 439 this_field%phase = -99d0 440 else 441 call errquit (pname//"invalid type: "//trim(type), 0, 0) 442 endif 443 444 445C XXX REFACTOR THIS 446 if (type .ne. "file") then 447 if (polarization.eq."x") then !not used yet, instead cases are done in excitation routine 448 this_field%theta = 90d0 449 this_field%phi = 0d0 450 elseif (polarization.eq."y") then ! "" 451 this_field%theta = 90d0 452 this_field%phi = 90d0 453 elseif (polarization.eq."z") then ! "" 454 this_field%theta = 0d0 455 this_field%phi = 0d0 456 elseif (polarization.eq."angle") then ! we actually use this 457 this_field%theta = theta 458 this_field%phi = phi 459 else 460 call errquit (pname//"invalid polarization: " 461 $ //trim(polarization), 0, 0) 462 endif 463 endif 464 465 call rt_tddft_field_rtdb_put (rtdb, nfields, this_field) 466 467 end subroutine 468 469 470C==================================================================== 471C 472C Generate entry name for field rtdb stuff (hack) 473C 474 subroutine rt_tddft_field_rtdb_entry_name (i, name) 475 implicit none 476 477#include "errquit.fh" 478#include "rtdb.fh" 479#include "mafdecls.fh" 480#include "stdio.fh" 481#include "rt_tddft.fh" 482 483 484C == Inputs == 485 integer, intent(in) :: i 486 487 488C == Outputs == 489 character(len=*), intent(out) :: name !was 17 490 491 492C == Parameters == 493 character(len=*), parameter :: pname = 494 $ "rt_tddft_field_rtdb_entry_name" 495 496 497C == Variables == 498 character*5 istring !note length 5 limit size of int 499 500 501 if ( (i .gt. rt_max_fields).or.(i .lt. 1) ) 502 $ call errquit(pname//"i must be between 1, rt_max_fields",0,0) 503 504 if (rt_max_fields .gt. 999) call errquit(pname// 505 $ "rt_max_fields too large; fix formatting", 0, 0) 506 507 write (istring, "(i0.5)") i 508 509 name = "rt_tddft:field_"//trim(istring)//"_" 510 511 end subroutine 512 513 514C==================================================================== 515C 516C Load field into rtbd. This is an ugly hack, but it's easier than 517C adding a custom struct to the rtdb routines. 518C 519 subroutine rt_tddft_field_rtdb_put (rtdb, i, field) 520 implicit none 521 522#include "rt_tddft.fh" 523#include "errquit.fh" 524#include "rtdb.fh" 525#include "mafdecls.fh" 526#include "stdio.fh" 527 528 529C == Inputs == 530 integer, intent(in) :: rtdb 531 integer, intent(in) :: i !index for the field 532 type(rt_field_t), intent(in) :: field 533 534 535C == Parameters == 536 character(len=*), parameter :: pname = "rt_tddft_field_rtdb_put: " 537 538 539C == Variables == 540 character*32 basename 541 character*32 entry_name 542 543 if ( (i .gt. rt_max_fields).or.(i .lt. 1) ) 544 $ call errquit(pname//"i must be between 1, rt_max_fields",0,0) 545 546 call rt_tddft_field_rtdb_entry_name (i, basename) 547 548 entry_name = trim(basename) // "name" 549 if (.not.rtdb_cput(rtdb,entry_name,1,field%name)) 550 $ call errquit(pname//'Write failed to name rtdb', 551 $ 0,RTDB_ERR) 552 553 entry_name = trim(basename) // "filename" 554 if (.not.rtdb_cput(rtdb,entry_name,1,field%filename)) 555 $ call errquit(pname//'Write failed to filename rtdb', 556 $ 0,RTDB_ERR) 557 558 entry_name = trim(basename) // "type" 559 if (.not.rtdb_cput(rtdb,entry_name,1,field%type)) 560 $ call errquit(pname//'Write failed to type rtdb', 561 $ 0,RTDB_ERR) 562 563 entry_name = trim(basename) // "polarization" 564 if (.not.rtdb_cput(rtdb,entry_name,1,field%polarization)) 565 $ call errquit(pname//'Write failed to polarization rtdb', 566 $ 0,RTDB_ERR) 567 568 entry_name = trim(basename) // "spin" 569 if (.not.rtdb_cput(rtdb,entry_name,1,field%spin)) 570 $ call errquit(pname//'Write failed to spin rtdb', 571 $ 0,RTDB_ERR) 572 573 entry_name = trim(basename) // "max" 574 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%max)) 575 $ call errquit(pname//'Write failed to max rtdb',0,RTDB_ERR) 576 577 entry_name = trim(basename) // "frequency" 578 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%frequency)) 579 $ call errquit(pname//'Write failed to frequency rtdb', 580 $ 0,RTDB_ERR) 581 582 entry_name = trim(basename) // "width" 583 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%width)) 584 $ call errquit(pname//'Write failed to width rtdb', 585 $ 0,RTDB_ERR) 586 587 entry_name = trim(basename) // "center" 588 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%center)) 589 $ call errquit(pname//'Write failed to center rtdb', 590 $ 0,RTDB_ERR) 591 592 entry_name = trim(basename) // "start" 593 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%start)) 594 $ call errquit(pname//'Write failed to start rtdb', 595 $ 0,RTDB_ERR) 596 597 entry_name = trim(basename) // "phase" 598 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%phase)) 599 $ call errquit(pname//'Write failed to phase rtdb', 600 $ 0,RTDB_ERR) 601 602 entry_name = trim(basename) // "theta" 603 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%theta)) 604 $ call errquit(pname//'Write failed to theta rtdb', 605 $ 0,RTDB_ERR) 606 607 entry_name = trim(basename) // "phi" 608 if (.not.rtdb_put(rtdb,entry_name,mt_dbl,1,field%phi)) 609 $ call errquit(pname//'Write failed to phi rtdb', 610 $ 0,RTDB_ERR) 611 612 end subroutine 613 614 615 616C 617C Get field from rtdb and put into struct 618C 619 subroutine rt_tddft_field_rtdb_get (rtdb, i, field) 620 implicit none 621 622#include "rt_tddft.fh" 623#include "errquit.fh" 624#include "rtdb.fh" 625#include "mafdecls.fh" 626#include "stdio.fh" 627 628 629C == Inputs == 630 integer, intent(in) :: rtdb 631 integer, intent(in) :: i !index for the field 632 633 634C == Outputs == 635 type(rt_field_t), intent(out) :: field 636 637 638 639C == Parameters == 640 character(len=*), parameter :: pname = "rt_tddft_field_rtdb_get: " 641 642 643C == Variables == 644 character*32 basename 645 character*32 entry_name 646 647 648 if ( (i .gt. rt_max_fields).or.(i .lt. 1) ) 649 $ call errquit(pname//"i must be between 1, rt_max_fields",0,0) 650 651 call rt_tddft_field_rtdb_entry_name (i, basename) 652 653 654 entry_name = trim(basename) // "name" 655 if (.not.rtdb_cget(rtdb,entry_name,1,field%name)) 656 $ call errquit(pname//'Read failed for name rtdb', 657 $ 0,RTDB_ERR) 658 659 entry_name = trim(basename) // "filename" 660 if (.not.rtdb_cget(rtdb,entry_name,1,field%filename)) 661 $ call errquit(pname//'Read failed for name rtdb', 662 $ 0,RTDB_ERR) 663 664 entry_name = trim(basename) // "type" 665 if (.not.rtdb_cget(rtdb,entry_name,1,field%type)) 666 $ call errquit(pname//'Read failed for type rtdb', 667 $ 0,RTDB_ERR) 668 669 entry_name = trim(basename) // "polarization" 670 if (.not.rtdb_cget(rtdb,entry_name,1,field%polarization)) 671 $ call errquit(pname//'Read failed for polarization rtdb', 672 $ 0,RTDB_ERR) 673 674 entry_name = trim(basename) // "spin" 675 if (.not.rtdb_cget(rtdb,entry_name,1,field%spin)) 676 $ call errquit(pname//'Read failed for spin rtdb', 677 $ 0,RTDB_ERR) 678 679 entry_name = trim(basename) // "max" 680 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%max)) 681 $ call errquit(pname//'Read failed for max rtdb',0,RTDB_ERR) 682 683 entry_name = trim(basename) // "frequency" 684 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%frequency)) 685 $ call errquit(pname//'Read failed for frequency rtdb', 686 $ 0,RTDB_ERR) 687 688 entry_name = trim(basename) // "width" 689 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%width)) 690 $ call errquit(pname//'Read failed for width rtdb', 691 $ 0,RTDB_ERR) 692 693 entry_name = trim(basename) // "center" 694 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%center)) 695 $ call errquit(pname//'Read failed for center rtdb', 696 $ 0,RTDB_ERR) 697 698 entry_name = trim(basename) // "start" 699 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%start)) 700 $ call errquit(pname//'Read failed for start rtdb', 701 $ 0,RTDB_ERR) 702 703 entry_name = trim(basename) // "phase" 704 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%phase)) 705 $ call errquit(pname//'Read failed for phase rtdb', 706 $ 0,RTDB_ERR) 707 708 entry_name = trim(basename) // "theta" 709 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%theta)) 710 $ call errquit(pname//'Read failed for theta rtdb', 711 $ 0,RTDB_ERR) 712 713 entry_name = trim(basename) // "phi" 714 if (.not.rtdb_get(rtdb,entry_name,mt_dbl,1,field%phi)) 715 $ call errquit(pname//'Read failed for phi rtdb', 716 $ 0,RTDB_ERR) 717 718 end subroutine 719c $Id$ 720