Added support for directly-interpolated waveforms

Added support for VP_SHAPE='alt spline', which directly interpolates waveforms
for the vector-potential. This form is intended as a replacement for the more
restrictive VP_SHAPE='table', and should be particularly useful for calculating
RABBITT and streaking spectra.

Also added new config files for gfortran 12 and Intel OneAPI compilers, and
made some cosmetic changes to the makefiles.

Author: Serguei Patchkovskii

Diff for: .gitignore

@@ -3,6 +3,7 @@
 *.x
 *.il
 *.a
+*.swp
 *__genmod.f90 
 *.optrpt
 *.dbg

Diff for: Makefile

@@ -1,3 +1,5 @@
+.PHONY: goal clear
+
 goal: makefile.dep
 	  $(MAKE) spherical_tdse.x
 	  # make build_pes.x
@@ -22,6 +24,9 @@ ACT2 = -e 's/^!\*nm/    /' # Disable MPI statements
 # include configs/babel-ifort18_opt.mak
 # include configs/zen-gfortran-7_opt.mak
 # include configs/zen-gfortran-11_opt.mak
+# include configs/zen-gfortran-12_opt.mak
+# include configs/zen-oneapi_opt.mak
+# include configs/zen-oneapi_opt_mpi.mak
 # include configs/zen-aocc-1.1_opt.mak      # VERY SLOW CODE. DO NOT USE.
 # include configs/oink-gfortran_opt.mak
 # include configs/macos_m1-gfortran_opt.mak

Diff for: configs/zen-gfortran-12_opt.mak

@@ -0,0 +1,17 @@
+BUILD_ID :="Optimized gfortran-12 (ZEN/EPYC), built on $(shell hostname) at $(shell date)"
+ACT = sed -e 's/^!\*qd/    /' # Enable quad-math statements
+# WARNING: requesting -mtune=znver2 leads to incorrect code
+# WARNING: requesting -fexternal-blas leads to unreasonable stack requirements
+F90 = gfortran-12 -I. \
+      -flto -O3 -fprotect-parens -march=native -mtune=native -fopenmp \
+      -ffast-math -fcx-fortran-rules -mrecip \
+      -floop-block \
+      -fno-realloc-lhs -fbacktrace -g \
+      -static-libgfortran \
+      -cpp -D__BUILD_ID__='$(BUILD_ID)' -ffree-line-length-none
+#     -fexternal-blas -fblas-matmul-limit=50 # Causes memory errors
+#     -floop-block # Requires Graphite isl, not configured on many distributions
+#     -ffast-math -fcx-fortran-rules -mrecip 
+F90L = $(F90) 
+LAPACK = -lopenblas_openmp # -lquadlapack_gfortran
+LIBEXTRA = -B/usr/share/libhugetlbfs -lhugetlbfs -Wl,--hugetlbfs-align

Diff for: configs/zen-oneapi_opt.mak

@@ -0,0 +1,17 @@
+BUILD_ID :="Optimized Intel oneAPI, built on $(shell hostname) at $(shell date)"
+ACT = sed -e 's/^!\*qd/    /' # Enable quad-math statements
+#ACT2 = -e 's/^!\*mp/    /' # Enable MPI statements
+F90 = ifort \
+            -qopenmp -qmkl=sequential -align all -align array256byte -pad \
+            -warn -assume buffered_io \
+            -O3 -ipo8 -no-prec-div -xAVX -axAVX2,MIC-AVX512 -complex_limited_range -fp-model fast=1 -ftz -assume protect_parens -heap-arrays 32 \
+            -qopt-prefetch=0 -qopt-matmul -qopt-subscript-in-range -qopt-dynamic-align \
+            -qopt-mem-layout-trans=3 -qopt-multi-version-aggressive \
+            -debug full -debug extended -traceback -qopt-report=5 \
+            -cpp -D__BUILD_ID__='$(BUILD_ID)'
+F90L = $(F90) 
+LAPACK = # Lapack is in Intel MKL, 
+LIBEXTRA = \
+           -static-intel -qopenmp-link=static \
+           -lhugetlbfs -Wl,-z,common-page-size=2097152 -Wl,-z,max-page-size=2097152 \
+           -Wl,--start-group $(MKLROOT)/lib/intel64/libmkl_intel_lp64.a $(MKLROOT)/lib/intel64/libmkl_core.a $(MKLROOT)/lib/intel64/libmkl_sequential.a -Wl,--end-group -lpthread -lm 

Diff for: configs/zen-oneapi_opt_mpi.mak

@@ -0,0 +1,17 @@
+BUILD_ID :="Optimized Intel oneAPI (MPI), built on $(shell hostname) at $(shell date)"
+ACT = sed -e 's/^!\*qd/    /' # Enable quad-math statements
+ACT2 = -e 's/^!\*mp/    /' # Enable MPI statements
+F90 = mpiifort \
+            -qopenmp -qmkl=sequential -align all -align array256byte -pad \
+            -warn -assume buffered_io \
+            -O3 -ipo8 -no-prec-div -xAVX -axAVX2,MIC-AVX512 -complex_limited_range -fp-model fast=1 -ftz -assume protect_parens -heap-arrays 32 \
+            -qopt-prefetch=0 -qopt-matmul -qopt-subscript-in-range -qopt-dynamic-align \
+            -qopt-mem-layout-trans=3 -qopt-multi-version-aggressive \
+            -debug full -debug extended -traceback -qopt-report=5 \
+            -cpp -D__BUILD_ID__='$(BUILD_ID)'
+F90L = $(F90) 
+LAPACK = # Lapack is in Intel MKL, 
+LIBEXTRA = \
+           -static-intel -static_mpi -qopenmp-link=static \
+           -lhugetlbfs -Wl,-z,common-page-size=2097152 -Wl,-z,max-page-size=2097152 \
+           -Wl,--start-group $(MKLROOT)/lib/intel64/libmkl_intel_lp64.a $(MKLROOT)/lib/intel64/libmkl_core.a $(MKLROOT)/lib/intel64/libmkl_sequential.a -Wl,--end-group -lpthread -lm 

Diff for: doc/CHANGES.txt

@@ -1,9 +1,13 @@
-Last updated: 2022 Jul 19
+Last updated: 2022 September 26
 ------------
 
 Change log
 ----------
 
+2022 Sep 26
+
+ - Added vp_shape="alt spline".
+
 2022 Jul 19
 
  - Added vp_shape="zz Sin2".

Diff for: doc/INPUT.txt

@@ -1,4 +1,4 @@
-Last updated: 2021 Aug 10
+Last updated: 2022 September 26
 ------------
 
   All input to scid-tdse is provided through the sph_tdse namelist, with one
@@ -990,10 +990,23 @@ Input keywords
                       description of the format. VP_SCALE has no effect.
                       WARNING: VP_SHAPE="table" is powerful, but makes it easy to 
                       WARNING: request a non-sensical calculation. Here be dragons.
+                      For most cases, VP_SHAPE="alt spline" is probably a better
+                      choice, as it is less dependent on the parameters of the 
+                      calculation.
       'spline'      - The vector potential is defined by cubic splines, specified
                       in file given by VP_TABLE (see VP_TABLE for the description).
-                      The result of the spline interpolation will be multiplied by
+                      The vector-potential is given by envelope + phase. The 
+                      result of the spline interpolation will be multiplied by
                       VP_SCALE.
+      'alt spline'  - The vector-potential is defined by cubic splines, specified
+                      in VP_TABLE (see below for the description). The splines 
+                      interpolate the vector-potential directly, so that a much
+                      denser set of spline points will likely be needed than for
+                      'spline'. The interpolated vector-potential is scaled by
+                      VP_SCALE and is delayed in time by VP_PARAM(3). If two sets
+                      of splines are given in the VP_TABLE, the second set is
+                      scaled and delayed by VP_SCALE_X and VP_PARAM_X(3), 
+                      respectively.
     "Gaussian" pulses above have envelopes modified at the wings to have 
     finite duration. See code in GaussianVP() in vectorpotential_tools.f90 
     for further details. The parameters are (add _X for the secondary component):
@@ -1016,11 +1029,14 @@ Input keywords
       VP_PARAM(4)  - Full width at zero [atomic units]
     "Flat-Sin2" pulses take the same parameters as "Sin2", plus:
       VP_PARAM(11) - Duration of the raising/falling edges [atomic units]
+    "alt spline" pulses depend on the following parameters:
+      VP_SCALE     - overall vector potential magnitude
+      VP_PARAM(3)  - Pulse time offset [atomic units]
     Also see FIELD_PREVIEW.
 
   VP_TABLE='vp.table' (vectorpotential_tools.f90)
-    Name of the file containing defining the vector-potential. Relevant 
-    if VP_SHAPE='table' and VP_SHARE='spline'. Also pay attention to the
+    Name of the file containing defining the vector-potential. Relevant if 
+    VP_SHAPE='table', 'spline', and 'alt spline'. Also pay attention to the
     value of ROTATION_MODE, as discussed below.
 
     For VP_SHAPE='table', the file must contain the contents of the 
@@ -1053,19 +1069,30 @@ Input keywords
     that all polar angles are zero and set ROTATION_MODE='none'. This also 
     disables vector-potential unwrapping.
 
-    For VP_SHAPE='spline', the input file is expected to defined six splines,
+    For VP_SHAPE='spline', the input file is expected to contain six splines,
     two for each Cartesian direction. For each pair, the first spline defines
     the Cartesian amplitude SA_{x/y/z}(t) of the vector-potential. The second 
     spline gives the corresponding phase SP_{x/y/z}(t), so that the vector 
     potential along that direction is given by:
 
         A_{x/y/z}(t) = SA_{x/y/z}(t) * cos(SP_{x/y/z}(t))
 
-    The resulting vector-potential will undergo field-unwrapping (unless
-    disabled by FIELD_UNWRAP=.FALSE.), and will be assumed to have general
-    polarization (unless linear polarization along the Z axis is explicitly
-    requested, by setting ROTATION_MODE='none'). It is expected that the
-    vector-potential is physical; no attempt to enforce this is made.
+    For VP_SHAPE='alt spline', the input file is expected to contain six
+    splines. The first three splines (one for each Cartesian direction) are
+    delayed in time by VP_PARAM(3) and scaled by VP_SCALE. The following 
+    three splines are delayed by VP_PARAM_X(3) and scaled by VP_SCALE_X, so
+    that the vector potential along Cartesian directions is given by:
+
+        A_x(t) = VP_SCALE*SP_1(t-VP_PARAM(3)) + VP_SCALE_X*SP_4(t-VP_PARAM_X(3))
+        A_y(t) = VP_SCALE*SP_2(t-VP_PARAM(3)) + VP_SCALE_X*SP_5(t-VP_PARAM_X(3))
+        A_z(t) = VP_SCALE*SP_3(t-VP_PARAM(3)) + VP_SCALE_X*SP_6(t-VP_PARAM_X(3))
+
+    For both VP_SHAPE='spline' and 'alt spline', the resulting vector-potential
+    will undergo field-unwrapping (unless disabled by FIELD_UNWRAP=.FALSE.), 
+    and will be assumed to have general polarization (unless linear polarization
+    along the Z axis is explicitly requested, by setting ROTATION_MODE='none'). 
+    It is expected that the vector-potential is physical; no attempt to enforce 
+    this is made.
 
     For each spline, we expect:
 

Diff for: examples/Makefile

@@ -0,0 +1,19 @@
+.PHONY: clean single cheap medium all
+
+all:
+	./run-all.sh all
+
+clean:
+	./clean.sh
+
+single:
+	./run-all.sh single hydrogen_1S_2P0_uniform.inp
+
+cheap:
+	./run-all.sh cheap
+
+medium:
+	./run-all.sh medium
+
+all:
+	./run-all.sh all

Diff for: math.f90

@@ -36,7 +36,7 @@ module math
   public MathInterpolate
   public rcsid_math
 !
-  character(len=clen), save :: rcsid_math = "$Id: math.f90,v 1.12 2021/04/26 15:44:44 ps Exp ps $"
+  character(len=clen), save :: rcsid_math = "$Id: math.f90,v 1.13 2022/02/19 18:54:26 ps Exp $"
 !
 !
   integer(ik), parameter       :: factorial_slack = 5         ! Extra factorials to produce while filling the cache

Diff for: vectorpotential_tools.f90

@@ -1,6 +1,6 @@
 !
 !   SCID-TDSE: Simple 1-electron atomic TDSE solver
-!   Copyright (C) 2015-2021 Serguei Patchkovskii, [email protected]
+!   Copyright (C) 2015-2022 Serguei Patchkovskii, [email protected]
 !
 !   This program is free software: you can redistribute it and/or modify
 !   it under the terms of the GNU General Public License as published by
@@ -59,7 +59,7 @@ module vectorpotential_tools
   public vp_apot
   public rcsid_vectorpotential_tools
   !
-  character(len=clen), save :: rcsid_vectorpotential_tools = "$Id: vectorpotential_tools.f90,v 1.19 2022/07/19 10:03:20 ps Exp $"
+  character(len=clen), save :: rcsid_vectorpotential_tools = "$Id: vectorpotential_tools.f90,v 1.21 2022/09/26 17:19:41 ps Exp ps $"
   !
   !  List of vector-potential names. These are used in a couple different places; 
   !  I would prefer any typos to cause a compile-time error!
@@ -78,6 +78,7 @@ module vectorpotential_tools
   character(len=*), parameter :: vpn_xySin2     = 'xy Sin2'
   character(len=*), parameter :: vpn_zxCW       = 'zx CW'
   character(len=*), parameter :: vpn_spline     = 'spline'
+  character(len=*), parameter :: vpn_altspline  = 'alt spline'
   !
   real(xk), save            :: vp_scale   = 0.00_xk     ! Overall magnitude of the vector potential
   character(len=clen), save :: vp_shape   = ' '         ! Shape of the vector-potential. See vp_apot()
@@ -105,6 +106,11 @@ module vectorpotential_tools
                                                         ! is described in doc/INPUT.txt, under the vp_table keyword.
                                                         ! vp_scale is influential in this case; don't forget to set it.
                                                         !
+                                                        ! For vp_shape='alt spline', the file should contain six individual
+                                                        ! splines, which are interpreted as two vector-potentials shifted 
+                                                        ! in time by VP_PARAM(3)/VP_PARAM_X(3), and scaled by 
+                                                        ! VP_SCALE/VP_SCALE_X. See doc/INPUT.txt for details
+                                                        !
                                                         ! In either case, the data should be supplied in Fortran 
                                                         ! free-form input format, with no fixed columns or formats.
   !
@@ -141,12 +147,16 @@ module vectorpotential_tools
   real(xk), save            :: flat_raise   = 0._xk     ! vp_param(11) : Duration of the raising and falling edges, in 
                                                         !                atomic units of time
   !
-  !  Data for vp_shape='spline', loaded from the file specified by vp_table
+  !  Data for vp_shape='spline' or 'alt spline', loaded from the file specified by vp_table
   !
-  type(cs_data), save       :: splines(6)               ! We define six splines in total, two for each Cartesian direction:
+  type(cs_data), save       :: splines(6)               ! vp_shape='spline'
+                                                        ! We define six splines in total, two for each Cartesian direction:
                                                         ! 1,2: Envelope and phase for A_x
                                                         ! 3,4: Ditto, for A_y
                                                         ! 5,6: Ditto, for A_z
+                                                        ! vp_shape='alt spline'
+                                                        ! 1,2,3: X/Y/Z components of the vector-potential 1
+                                                        ! 4,5,6: X/Y/Z components of the vector-potential 2
   !
   contains
   !
@@ -252,6 +262,8 @@ function vp_apot(t,theta,phi) result(apot)
         ph = atan2(ay,ax)
       case (vpn_spline)
         call vp_spline(t,apot,th,ph)
+      case (vpn_altspline)
+        call vp_altspline(t,apot,th,ph)
     end select
     if (present(theta)) theta = th
     if (present(phi)  ) phi   = ph
@@ -301,6 +313,8 @@ subroutine initialize
           omega = vp_param(1)
         case (vpn_spline)
           call init_splines
+        case (vpn_altspline)
+          call init_altsplines
       end select
       first = .false.
       !$omp flush(first)
@@ -523,9 +537,6 @@ subroutine vp_spline(t,vpot,th,ph)
   end subroutine vp_spline
   !
   subroutine init_splines
-    integer(ik)                  :: iu, ios, ispl, npts
-    real(xk), allocatable        :: xy(:,:)
-    character(len=clen)          :: msg
     character(len=10), parameter :: names(6) = (/ 'X envelope', 'X phase   ', &
                                                   'Y envelope', 'Y phase   ', &
                                                   'Z envelope', 'Z phase   ' /)
@@ -534,6 +545,56 @@ subroutine init_splines
     write (out,"(t5,'Loading splines from ',a)") trim(vp_table)
     write (out,"(t5,'Scaling factor = ',g0.8)") vp_scale
     !
+    call init_splines_common(names)
+  end subroutine init_splines
+  !
+  subroutine vp_altspline(t,vpot,th,ph)
+    real(xk), intent(in)  :: t            ! Time where we need the vector-potential
+    real(xk), intent(out) :: vpot, th, ph ! Vector-potential, spherical coordinates
+    !
+    real(xk) :: vp_x, vp_y, vp_z ! Cartesian components of the vector-potential
+    !
+    vp_x =        vp_scale   * cs_evaluate(splines(1),t-vp_param  (3))
+    vp_y =        vp_scale   * cs_evaluate(splines(2),t-vp_param  (3))
+    vp_z =        vp_scale   * cs_evaluate(splines(3),t-vp_param  (3))
+    vp_x = vp_x + vp_scale_x * cs_evaluate(splines(4),t-vp_param_x(3))
+    vp_y = vp_y + vp_scale_x * cs_evaluate(splines(5),t-vp_param_x(3))
+    vp_z = vp_z + vp_scale_x * cs_evaluate(splines(6),t-vp_param_x(3))
+    !
+    !  Special case: if the X and Y fields vanish identically, force theta and phi to be zero
+    !  In all other cases, the unwrap routine should take of keeping the vector-potential
+    !  as smooth as possible.
+    !
+    if (vp_x==0._xk .and. vp_y==0._xk) then
+      vpot = vp_z ; th = 0._xk ; ph = 0._xk
+    else
+      vpot = sqrt(vp_x**2+vp_y**2+vp_z**2)
+      th   = acos(vp_z/vpot)
+      ph   = atan2(vp_y,vp_x)
+    end if
+  end subroutine vp_altspline
+  !
+  subroutine init_altsplines
+    character(len=2), parameter :: xyz(6) = (/ 'X1', 'Y1', 'Z1', 'X2', 'Y2', 'Z2' /)
+    !
+    write (out,"(/'Using cubic-spline interpolation to define the vector-potential')")
+    write (out,"(t5,'Loading splines from ',a)") trim(vp_table)
+    write (out,"(t5,'Scaling factor [Field 1] = ',g24.12)") vp_scale
+    write (out,"(t5,'    Time delay [Field 1] = ',g24.12)") vp_param(3)
+    write (out,"(t5,'Scaling factor [Field 2] = ',g24.12)") vp_scale_x
+    write (out,"(t5,'    Time delay [Field 2] = ',g24.12)") vp_param_x(3)
+    write (out,"()")
+    !
+    call init_splines_common(xyz)
+  end subroutine init_altsplines
+
+  subroutine init_splines_common(codes)
+    character(len=*), intent(in) :: codes(:)
+    !
+    integer(ik)                 :: iu, ios, ispl, npts
+    real(xk), allocatable       :: xy(:,:)
+    character(len=clen)         :: msg
+    !
     error_block: do
       msg = "opening " // trim(vp_table)
       open(newunit=iu,form='formatted',recl=256,action='read',position='rewind',status='old', &
@@ -548,7 +609,7 @@ subroutine init_splines
         allocate (xy(2,npts),stat=ios)
         if (ios/=0) exit error_block
         !
-        write (out,"(t5,'spline ',i0,' defining ',a,' contains ',i0,' points.')") ispl, trim(names(ispl)), npts
+        write (out,"(t5,'spline ',i0,' defining component ',a,' contains ',i0,' points.')") ispl, trim(codes(ispl)), npts
         !
         if (npts>0) then
           write (msg,"('reading data points for spline #',i0,' npts = ',i0)") ispl, npts
@@ -568,7 +629,7 @@ subroutine init_splines
       return
     end do error_block
     !
-    write (out,"('init_splines: Error ',a,'. ios = ',i0)") trim(msg), ios
-    stop 'vectorpotential_tools%init_splines - error initializing splines'
-  end subroutine init_splines
+    write (out,"('init_splines_common: Error ',a,'. ios = ',i0)") trim(msg), ios
+    stop 'vectorpotential_tools%init_splines_common - error initializing splines'
+  end subroutine init_splines_common
 end module vectorpotential_tools