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