src/.ycm_extra_conf.py

@@ -36,6 +36,11 @@
     '-lmpi',
     '-lmpi_cxx',
 
+    '-I{}/gsl/include'.format(environ['QPROP_DEP_DIR']),
+    '-L{}/gsl/lib'.format(environ['QPROP_DEP_DIR']),
+    '-lgsl',
+    '-lgslcblas',
+
     # QPROP flags
     '-L{}/lib/x86_64'.format(environ['QPROP_HOME']),
     '-I{}/src/base'.format(environ['QPROP_HOME']),

src/base/grid.cc

@@ -1 +1,9 @@
 #include <grid.h>
+
+long num_of_basis(long qprop_dim, long ell_grid_size) {
+  switch (qprop_dim) {
+    case 34: return ell_grid_size;
+    case 44: return ell_grid_size * ell_grid_size;
+    default: return -1;
+  }
+}

src/base/grid.h

@@ -202,4 +202,9 @@ class grid {
   long offse_x, offse_y, offse_z;
 };
 
+
+
+long num_of_basis(long qprop_dim, long ell_grid_size);
+
+
 #endif 

src/base/potentials.hh

@@ -14,30 +14,41 @@ double always_zero_imag(long x, long y, long z, double t, grid g);
 
 // The vector potential with sine squared envelope
 class vecpot {
-  double omega;
-  double n_c;
-  double E_0;
-  double phi_cep;
-  double duration;
-  double ww;
+
+  double omega;  // angular freqency
+  double n_c;  // number of cycles
+  double E_0;  // max electric field amplitude
+  double phi_cep;  // carrier-envelop-phase (CEP)
+  double duration;  // duration of pulse
+  double ww;  // angular frequency of the envelope
   double start_time, end_time;
+
 public:
+
   vecpot() { duration = -1; };
-  vecpot(double om, double n_cyc, double E_max, double cep, double start_time_in = 0.0) : omega(om), n_c(n_cyc), E_0(E_max), phi_cep(cep), start_time(start_time_in) {
-    duration=n_c*2*M_PI/omega;
+
+  vecpot(
+      double om, double n_cyc, double E_max, double cep, 
+      double start_time_in = 0.0) : 
+    omega(om), n_c(n_cyc), E_0(E_max), phi_cep(cep), start_time(start_time_in) 
+  {
+    duration = n_c * 2*M_PI/omega;
     // angular frequency of the envelope
     ww=omega/(2.0*n_c);
     end_time = start_time + duration;
   };
+
   double operator()(double time, int me) const {
     // switch of the field for propagation after the pulse
-    if ((time>0.0) && (time<duration)) {
-      return E_0/omega*pow2(sin(ww*time))*sin(omega*time+phi_cep); // here's the shape of the laser pulse
+    if ((time>start_time) && (time<end_time)) {
+      // here's the shape of the laser pulse
+      return E_0/omega*pow2(sin(ww*(time-start_time)))*sin(omega*(time-start_time)+phi_cep); 
     }
     else {
-      return 0;
+      return 0.0;
     };
   };
+
   double get_omega() { return omega; }
   double get_E0() { return E_0; }
   double get_num_cycles() { return n_c; }
@@ -50,7 +61,7 @@ public:
   double get_Up() {
     return (E_0 * E_0)/4.0/(omega*omega);
   };
-  double integral(double Time) {
+  double integral(double absolute_time) {
     // grind(ratsimp(integrate(sin(pi*time/pulse_dur)**2*sin(omega*time+phi), time, 0, Time)));
     // const double omega = omega;
     const double omega_2 = pow2(omega);
@@ -63,7 +74,12 @@ public:
     const double cos_phi_2 = pow2(cos(phi_cep));
     const double phi = phi_cep;
     const double ampl = E_0/omega;
-    if (Time<pulse_dur) {
+
+    double Time = absolute_time - start_time;
+
+    if (Time < 0.0) { 
+      return 0.0; 
+    } else if (Time < duration) {
       return ampl*((omega_2*sin(phi)*pulse_dur_2-2.0*omega*sin(phi)*pi*pulse_dur)*sin(((omega*pulse_dur+2.0*pi)*Time+2.0*phi*pulse_dur)/pulse_dur)
 		   +(omega_2*cos(phi)*pulse_dur_2-2.0*omega*cos(phi)*pi*pulse_dur)*cos(((omega*pulse_dur+2.0*pi)*Time+2.0*phi*pulse_dur)/pulse_dur)
 		   +(omega_2*sin(phi)*pulse_dur_2+2.0*omega*sin(phi)*pi*pulse_dur)*sin(((omega*pulse_dur-2.0*pi)*Time+2.0*phi*pulse_dur)/pulse_dur)
@@ -78,8 +94,7 @@ public:
 		   +(8.0*cos(phi)*pi_2-2.0*omega_2*cos(phi)*pulse_dur_2)*cos(omega*Time)
 		   +(-8.0*sin(phi)*sin(2.0*phi)-8.0*cos(phi)*cos(2.0*phi)-8.0*cos(phi))*pi_2)
 	/((8.0*omega_3*sin_phi_2+8.0*omega_3*cos_phi_2)*pulse_dur_2+(-32.0*omega*sin_phi_2-32.0*omega*cos_phi_2)*pi_2);
-    }
-    else {
+    } else {
       return ampl*((omega_2*sin(phi)*pulse_dur_2-2.0*omega*sin(phi)*pi*pulse_dur)*sin(((omega*pulse_dur+2.0*pi)*pulse_dur+2.0*phi*pulse_dur)/pulse_dur)
 		   +(omega_2*cos(phi)*pulse_dur_2-2.0*omega*cos(phi)*pi*pulse_dur)*cos(((omega*pulse_dur+2.0*pi)*pulse_dur+2.0*phi*pulse_dur)/pulse_dur)
 		   +(omega_2*sin(phi)*pulse_dur_2+2.0*omega*sin(phi)*pi*pulse_dur)*sin(((omega*pulse_dur-2.0*pi)*pulse_dur+2.0*phi*pulse_dur)/pulse_dur)
@@ -223,7 +238,12 @@ class superposed_vecpot {
     double min_start_time = vp.get_start_time();
     for (vecpot_index = 1; vecpot_index < num_of_vecpot; vecpot_index++) {
       vp = vecpot_arr[vecpot_index];
-      if (min_start_time > vp.get_start_time()) { min_start_time = vp.get_start_time(); }
+
+//      cout << "[ LOG ][0] vp.start_time = " << vp.get_start_time() << std::endl;
+
+      if (min_start_time > vp.get_start_time()) { 
+        min_start_time = vp.get_start_time(); 
+      }
     }
     return min_start_time;
   }
@@ -234,7 +254,9 @@ class superposed_vecpot {
     double max_end_time = vp.get_end_time();
     for (vecpot_index = 1; vecpot_index < num_of_vecpot; vecpot_index++) {
       vp = vecpot_arr[vecpot_index];
-      if (max_end_time < vp.get_end_time()) { max_end_time = vp.get_end_time(); }
+      if (max_end_time < vp.get_end_time()) { 
+        max_end_time = vp.get_end_time(); 
+      }
     }
     return max_end_time;
   }

src/base/vecpot.cc

@@ -6,7 +6,7 @@ using std::cout;
 
 bool is_valid_direction(char direction) {
     int char_index, num_of_matches = 0;
-    for (char_index=0; char_index<strlen(valid_directions); char_index++) {
+    for (char_index=0; char_index<int(strlen(valid_directions)); char_index++) {
 	if (direction == valid_directions[char_index]) { num_of_matches++; }
     }
     if ( num_of_matches == 1 ) { return true; } // there is a single matched valid direction
@@ -22,6 +22,7 @@ void print_vecpot_param(struct vecpot_param vparam) {
     cout << "E0: " << vparam.E0 << endl;
     cout << "num-cycles: " << vparam.num_cycles << endl;
     cout << "phase / pi " << vparam.phase_pi << endl;
+    cout << "start-time " << vparam.start_time << endl;
 }
 
 bool if_exist_get_vecpot_param ( parameterListe& para_prop, char direction, 
@@ -30,7 +31,9 @@ bool if_exist_get_vecpot_param ( parameterListe& para_prop, char direction,
   if ( ! is_valid_direction(direction) ) { 
 	  fprintf(stderr, "[ LOG ] The given direction('%c') is not valid.\n", direction);
 	  fprintf(stderr, "[ LOG ] Please choose direction among the following characters: \"%s\"\n", valid_directions);
-	  throw "[ERROR] Direction is not valid!"; }
+	  throw "[ERROR] Direction is not valid!"; 
+  }
+
   // construct suffix for each parameter name in a form: '-i-j'
   string suffix = string("-") + string(1, direction)
     + string("-") + to_string((long long int) vecpot_index); // 'long long int' is for working around ambiguous to_string method due to overloading
@@ -51,6 +54,13 @@ bool if_exist_get_vecpot_param ( parameterListe& para_prop, char direction,
   } catch (std::exception&) {
     vecpot_set_exist = false;
   }
+
+  try {
+    vp->start_time = para_prop.getDouble(string("start-time")+suffix);
+  } catch (std::exception&) {
+    cout << "[ LOG ] start-time set to default value.\n";
+  }
+
 
   if (vecpot_set_exist) { print_vecpot_param(*vp); }
 
@@ -92,6 +102,9 @@ int vecpot_param_with_original_param_34 ( parameterListe& para, struct vecpot_pa
 
   try { vp->phase_pi = para.getDouble("phase") / M_PI; }
   catch (std::exception&) { vp->phase_pi = 0.0; }  // default value
+
+  try { vp->start_time = para.getDouble("start-time"); }
+  catch (std::exception&) {} // The default value is defined in the vp object
 
   return 0;
 }
@@ -106,7 +119,8 @@ int construct_vecpot_with_original_param_34 ( parameterListe& para,
 
   vecpot_x = vecpot(vp.omega, 1.0, 0.0, 0.0);
   vecpot_y = vecpot(vp.omega, 1.0, 0.0, 0.0);
-  vecpot_z = vecpot(vp.omega, vp.num_cycles, vp.E0, vp.phase_pi * M_PI);
+  vecpot_z = vecpot(vp.omega, vp.num_cycles, vp.E0, vp.phase_pi * M_PI, 
+      vp.start_time);
 
   return 0;
 }
@@ -129,7 +143,10 @@ int if_exist_get_superposed_vecpot(parameterListe& para, char direction,
         << " and index " << vecpot_index_in_param_name << endl;
       return 1;
     }
-    vp_arr[vecpot_index] = vecpot(vparam.omega, vparam.num_cycles, vparam.E0, vparam.phase_pi * M_PI);
+    // Construct vecpot object
+    vp_arr[vecpot_index] = vecpot(
+        vparam.omega, vparam.num_cycles, vparam.E0, 
+        vparam.phase_pi * M_PI, vparam.start_time);
   }
   svp = superposed_vecpot(vp_arr, num_of_vecpot);
   return 0;
@@ -142,6 +159,7 @@ void print_vecpot(vecpot& vp, const char *tag) {
   cout << "[ LOG ] - E0 = " << vp.get_E0() << endl;
   cout << "[ LOG ] - num_cycles " << vp.get_num_cycles() << endl;
   cout << "[ LOG ] - phase " << vp.get_phase() << endl;
+  cout << "[ LOG ] - start_time " << vp.get_start_time() << endl;
   cout << "[ LOG ] ---------------------------------" << endl; 
 }
 
@@ -159,6 +177,7 @@ void print_superposed_vecpot(superposed_vecpot& svp, const char *tag) {
 
 int construct_superposed_vecpot_at_direction(parameterListe& para,
     char direction, superposed_vecpot& svp) {
+
   int number_of_vecpot_param_set = get_number_of_vecpot_param_set(para, direction);
   int return_code = if_exist_get_superposed_vecpot(para, 
       direction, svp, number_of_vecpot_param_set);
@@ -172,10 +191,10 @@ int construct_superposed_vecpot_at_direction(parameterListe& para,
 }
 
 
-//int construct_vecpot ( long dim, parameterListe& para_prop, 
-//    vecpot& vecpot_x, vecpot& vecpot_y, vecpot& vecpot_z ) {
 int construct_vecpot ( long dim, parameterListe& para_prop, 
-    superposed_vecpot& vecpot_x, superposed_vecpot& vecpot_y, superposed_vecpot& vecpot_z ) {
+    superposed_vecpot& vecpot_x, superposed_vecpot& vecpot_y, 
+    superposed_vecpot& vecpot_z ) 
+{
 
   if (dim==34) {
     //char direction = 'z';
@@ -237,5 +256,7 @@ int construct_vecpot ( long dim, parameterListe& para_prop,
     std::cerr << "[ERROR] Unexpected propagation mode (dimension): "
       << dim << endl;
   }
+
+  return 0;
 }
 

src/base/vecpot.hh

@@ -12,7 +12,16 @@
 
 
 struct vecpot_param {
-  double omega, E0, num_cycles, phase_pi;
+  double omega, E0, num_cycles, phase_pi, start_time;
+
+  vecpot_param(): start_time(0.0) {};
+
+  vecpot_param(
+      double omega_in, double E0_in, double num_cycles_in, double phase_pi_in, 
+      double start_time_in=0.0):
+    omega(omega_in), num_cycles(num_cycles_in), phase_pi(phase_pi_in),
+    start_time(start_time_in) {};
+
 };
 
 const char valid_directions[] = "xyz";

src/example/ati-tsurff-winop/initial.param

@@ -0,0 +1,8 @@
+nuclear-charge double 1.0
+pot-cutoff double 25.0
+delta-r double 0.2
+radial-grid-size double 100.0
+ell-grid-size long 1
+qprop-dim long 34
+initial-m long 0
+initial-ell long 0

src/example/ati-tsurff-winop/propagate.param

@@ -0,0 +1,7 @@
+imag-width double 150.0
+ell-grid-size long 15
+max-electric-field double 0.02387
+omega double 0.085
+num-cycles double 20.0
+# time step for propagation; delta-r/4 is a sensible choice
+delta-t double 0.05

src/example/ati-tsurff-winop/run.sh

@@ -0,0 +1,21 @@
+#!/bin/bash
+
+
+## Step 1
+$QPROP_HOME/bin/imag-prop
+$QPROP_HOME/bin/real-prop
+
+
+## Step 2 (select either one or both)
+
+# for tsurff
+$QPROP_HOME/bin/ppp
+$QPROP_HOME/bin/eval-tsurff
+cp tsurff-polar0.dat tsurff-polar.dat
+cp tsurff-partial0.dat tsurff-partial.dat
+
+# for winop
+$QPROP_HOME/bin/winop
+cp spectrum_0.dat spectrum.dat
+cp spectrum_polar0.dat spectrum_polar.dat
+

src/example/ati-tsurff-winop/tsurff.param

@@ -0,0 +1,21 @@
+# distance to the t-SURFF boundary
+R-tsurff double 100.0
+# This determines the duration of the simulation (slowest electron to reach the t-SURFF boundary).
+p-min-tsurff double 0.125
+# largest k value in the calculated spectrum
+k-max-surff double 1.2
+# number of k values for the spectrum
+num-k-surff long 101
+# number of angles theta (\theta \in [0:\pi])
+num-theta-surff long 51
+# number of angles phi (\phi \in [0:2\pi])
+num-phi-surff long 1
+# delta-k-scheme=1: equidistant k grid discretization w.r.t momentum; delta-k-scheme=2: equidistant k grid discretization w.r.t energy
+delta-k-scheme long 2
+# how many time steps are processed at a time during evaluation of the spectrum
+cache-size-t long 512
+# expansion-scheme=2: true expansion of the spectrum in spherical harmonics
+# expansion-scheme=1: use for small number of angles and large number of ells (no partial spectra are produced)
+expansion-scheme long 2
+
+use-ppp bool 1

src/example/ati-tsurff-winop/winop.param

@@ -0,0 +1,11 @@
+energy-max double 0.72
+energy-min double 0.0
+num-energy long 101
+num-theta long 51
+num-phi long 1
+# improve the angular resolved spectrum by  propagating longer
+extra-time-winop double 0.0
+# improve the energy resolution by using a large radial grid for the application of the window operator
+winop-radial-grid-size long 50000
+
+

src/example/attoclock-tsurff-winop/initial.param

@@ -0,0 +1,8 @@
+nuclear-charge double 1.0
+pot-cutoff double 25.0
+delta-r double 0.2
+radial-grid-size double 100.0
+ell-grid-size long 1
+qprop-dim long 44
+initial-m long 0
+initial-ell long 0

src/example/attoclock-tsurff-winop/propagate.param

@@ -0,0 +1,17 @@
+imag-width double 150.0
+ell-grid-size long 25
+delta-t double 0.05
+
+max-electric-field-x-1 double 0.0533799
+omega-x-1 double 0.114
+num-cycles-x-1 double 2.0
+phase-pi-x-1 double 0.0
+
+max-electric-field-y-1 double 0.0533799
+omega-y-1 double 0.114
+num-cycles-y-1 double 2.0
+phase-pi-y-1 double 0.5
+
+#max-electric-field double 0.0533799
+#omega double 0.114
+#num-cycles double 2.0

src/example/attoclock-tsurff-winop/run.sh

@@ -0,0 +1,21 @@
+#!/bin/bash
+
+
+## Step 1
+$QPROP_HOME/bin/imag-prop
+$QPROP_HOME/bin/real-prop
+
+
+## Step 2 (select either one or both)
+
+# for tsurff
+$QPROP_HOME/bin/ppp
+$QPROP_HOME/bin/eval-tsurff
+cp tsurff-polar0.dat tsurff-polar.dat
+cp tsurff-partial0.dat tsurff-partial.dat
+
+# for winop
+$QPROP_HOME/bin/winop
+cp spectrum_0.dat spectrum.dat
+cp spectrum_polar0.dat spectrum_polar.dat
+