Added aimd_static examples. May not be complete...

spalato · spalato · commit 56d1d7a98d9e · 2022-11-02T12:46:36.000+01:00
diff --git a/examples/aimd_static/cfg_0_10.yaml b/examples/aimd_static/cfg_0_10.yaml
@@ -0,0 +1,38 @@
+e_frame: 1.75
+rootname: sys_0_10_xcorr_0_10
+sigma: 0.02
+states_fn: sys_0_10.txt
+t1: [0.0, 3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0, 27.0, 30.0, 33.0, 36.0, 39.0,
+  42.0, 45.0, 48.0, 51.0, 54.0, 57.0, 60.0, 63.0, 66.0, 69.0, 72.0, 75.0, 78.0, 81.0,
+  84.0, 87.0, 90.0, 93.0, 96.0, 99.0, 102.0, 105.0, 108.0, 111.0, 114.0, 117.0, 120.0,
+  123.0, 126.0, 129.0, 132.0, 135.0, 138.0, 141.0, 144.0, 147.0, 150.0, 153.0, 156.0,
+  159.0, 162.0, 165.0, 168.0, 171.0, 174.0, 177.0, 180.0, 183.0, 186.0, 189.0, 192.0,
+  195.0, 198.0, 201.0, 204.0, 207.0, 210.0, 213.0, 216.0, 219.0, 222.0, 225.0, 228.0,
+  231.0, 234.0, 237.0, 240.0, 243.0, 246.0, 249.0, 252.0, 255.0, 258.0, 261.0, 264.0,
+  267.0, 270.0, 273.0, 276.0, 279.0, 282.0, 285.0, 288.0, 291.0, 294.0, 297.0, 300.0,
+  303.0, 306.0, 309.0, 312.0, 315.0, 318.0, 321.0, 324.0, 327.0, 330.0, 333.0, 336.0,
+  339.0, 342.0, 345.0, 348.0, 351.0, 354.0, 357.0, 360.0, 363.0, 366.0, 369.0, 372.0,
+  375.0, 378.0, 381.0, 384.0, 387.0, 390.0, 393.0, 396.0, 399.0, 402.0]
+t2: [0.0, 3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0, 27.0, 30.0, 33.0, 36.0, 39.0,
+  42.0, 45.0, 48.0, 51.0, 54.0, 57.0, 60.0, 63.0, 66.0, 69.0, 72.0, 75.0, 78.0, 81.0,
+  84.0, 87.0, 90.0, 93.0, 96.0, 99.0, 102.0, 105.0, 108.0, 111.0, 114.0, 117.0, 120.0,
+  123.0, 126.0, 129.0, 132.0, 135.0, 138.0, 141.0, 144.0, 147.0, 150.0, 153.0, 156.0,
+  159.0, 162.0, 165.0, 168.0, 171.0, 174.0, 177.0, 180.0, 183.0, 186.0, 189.0, 192.0,
+  195.0, 198.0, 201.0, 204.0, 207.0, 210.0, 213.0, 216.0, 219.0, 222.0, 225.0, 228.0,
+  231.0, 234.0, 237.0, 240.0, 243.0, 246.0, 249.0, 252.0, 255.0, 258.0, 261.0, 264.0,
+  267.0, 270.0, 273.0, 276.0, 279.0, 282.0, 285.0, 288.0, 291.0, 294.0, 297.0, 300.0,
+  303.0, 306.0, 309.0, 312.0, 315.0, 318.0, 321.0, 324.0, 327.0, 330.0, 333.0, 336.0,
+  339.0, 342.0, 345.0, 348.0, 351.0, 354.0, 357.0, 360.0, 363.0, 366.0, 369.0, 372.0,
+  375.0, 378.0, 381.0, 384.0, 387.0, 390.0, 393.0, 396.0, 399.0]
+t3: [0.0, 3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0, 27.0, 30.0, 33.0, 36.0, 39.0,
+  42.0, 45.0, 48.0, 51.0, 54.0, 57.0, 60.0, 63.0, 66.0, 69.0, 72.0, 75.0, 78.0, 81.0,
+  84.0, 87.0, 90.0, 93.0, 96.0, 99.0, 102.0, 105.0, 108.0, 111.0, 114.0, 117.0, 120.0,
+  123.0, 126.0, 129.0, 132.0, 135.0, 138.0, 141.0, 144.0, 147.0, 150.0, 153.0, 156.0,
+  159.0, 162.0, 165.0, 168.0, 171.0, 174.0, 177.0, 180.0, 183.0, 186.0, 189.0, 192.0,
+  195.0, 198.0, 201.0, 204.0, 207.0, 210.0, 213.0, 216.0, 219.0, 222.0, 225.0, 228.0,
+  231.0, 234.0, 237.0, 240.0, 243.0, 246.0, 249.0, 252.0, 255.0, 258.0, 261.0, 264.0,
+  267.0, 270.0, 273.0, 276.0, 279.0, 282.0, 285.0, 288.0, 291.0, 294.0, 297.0, 300.0,
+  303.0, 306.0, 309.0, 312.0, 315.0, 318.0, 321.0, 324.0, 327.0, 330.0, 333.0, 336.0,
+  339.0, 342.0, 345.0, 348.0, 351.0, 354.0, 357.0, 360.0, 363.0, 366.0, 369.0, 372.0,
+  375.0, 378.0, 381.0, 384.0, 387.0, 390.0, 393.0, 396.0, 399.0, 402.0]
+xcorr_fn: xcorr_0_10.txt
diff --git a/examples/aimd_static/run_third_order.ps1 b/examples/aimd_static/run_third_order.ps1
@@ -0,0 +1,10 @@
+Set-PSDebug -Trace 1
+$cfg=$args[0]
+$root=$(Select-String -Path $cfg -Pattern '^rootname:\W*(.*)$').Matches.Groups[1].value
+julia .\third_order_xc.jl $cfg
+ls ${root}_r[nr].bin | %{julia .\spec_ft.jl $cfg $_ $_.Name.Replace("_r", "_s")}
+ls ${root}_s[nr].bin | %{python .\plot_spec.py $cfg $_}
+julia .\spec_abs.jl $cfg $root ${root}_sa.bin
+#rm *_[sr][rn].bin
+python .\plot_spec.py $cfg ${root}_sa.bin
+#python .\beatmap.py $cfg ${root}_sa.bin
diff --git a/examples/aimd_static/spec_abs.jl b/examples/aimd_static/spec_abs.jl
@@ -0,0 +1,28 @@
+# fourier transform coherences to 2d spectra
+# usage: julia spec_ft.jl <config> <out>
+
+import YAML
+
+function run(cfgf, root, outf)
+    cfg = open(YAML.load, cfgf)
+    sz = tuple(map(length, cfg["t$i"] for i=1:3)...)
+    inf_r = root*"_sr.bin"
+    inf_n = root*"_sn.bin"
+    in_r = open(inf_r, "r+")
+    in_n = open(inf_n, "r+")
+    s_r = Mmap.mmap(in_r, Array{Complex128, 3}, sz)
+    s_n = Mmap.mmap(in_n, Array{Complex128, 3}, sz)
+    out = open(outf, "w+")
+    try
+        spec_abs = Mmap.mmap(out, Array{Complex128, 3}, sz)
+        spec_abs[:] = copy(s_n)
+        spec_abs[2:end,:,:] += flipdim(s_r[2:end,:,:], 1)
+    finally
+        close(out)
+        close(in_r)
+        close(in_n)
+    end
+end
+
+cfgf, root, outf = ARGS[1:3]
+run(cfgf, root, outf)
diff --git a/examples/aimd_static/spec_ft.jl b/examples/aimd_static/spec_ft.jl
@@ -0,0 +1,24 @@
+# fourier transform coherences to 2d spectra
+# usage: julia spec_ft.jl <config> <in> <out>
+
+import YAML
+
+function run(cfgf, infn, outf)
+    cfg = open(YAML.load, cfgf)
+    sz = tuple(map(length, cfg["t$i"] for i=1:3)...)
+    inf = open(infn, "r+")
+    sig = Mmap.mmap(inf, Array{Complex128, 3}, sz)
+    sig[1,:,:] *= 0.5
+    sig[:,:,1] *= 0.5
+    out = open(outf, "w+")
+    spec2d = Mmap.mmap(out, Array{Complex128, 3}, sz)
+    fill!(spec2d, 0)
+    spec2d[:] = fftshift(ifft(sig, (1,3)), (1,3))
+    close(out)
+    sig[1,:,:] *= 2
+    sig[:,:,1] *= 2
+    close(inf)
+end
+
+cfgf, infn, outf = ARGS[1:3]
+run(cfgf, infn, outf)
diff --git a/examples/aimd_static/sys_0_10.txt b/examples/aimd_static/sys_0_10.txt
@@ -0,0 +1,10 @@
+1.561819839679357891e+00 5.421701363481752356e-01
+1.611593186372745734e+00 5.131993707826494333e-01
+1.760613226452905300e+00 2.833814072832257502e-01
+1.794191783567134602e+00 2.761258042416964154e-01
+1.820523647294588709e+00 2.376053379524454001e-01
+1.852621442885769687e+00 2.412396021944760716e-01
+1.880126653306613438e+00 2.926030893243291020e-01
+1.913341082164327744e+00 4.056453378656108000e-01
+1.935171943887774360e+00 4.839319692370925852e-01
+1.968658517034066469e+00 3.471332581373224313e-01
diff --git a/examples/aimd_static/third_order_xc.jl b/examples/aimd_static/third_order_xc.jl
@@ -0,0 +1,103 @@
+import YAML
+using Mbo
+using ProgressMeter
+
+tag(i) = @sprintf "X%02d" i
+function parse_corr(fn)
+    readdlm(fn)
+end
+
+function run(args)
+t0 = now()
+cfgf = args[1]
+info("Loading parameters from $(cfgf)")
+cfg = open(YAML.load, cfgf)
+states_fn = cfg["states_fn"]
+cf_fn = cfg["xcorr_fn"]
+
+t1 = convert.(Float64, cfg["t1"])
+t2 = convert.(Float64, cfg["t2"])
+t3 = convert.(Float64, cfg["t3"])
+sigma = cfg["sigma"]
+
+info("Loading energies from $states_fn")
+states = readdlm(states_fn)
+energies = states[:,1]
+tdm = states[:,2]
+info("States count: $(size(states)[1])")
+info("Loading correlation matrix from $cf_fn")
+corrm = parse_corr(cf_fn)
+info("Correlation count: $(length(corrm))")
+
+# build system
+info("Building system...")
+t0_sys = now()
+s = System("G")
+frame = ev2angphz(cfg["e_frame"])
+lut_grid = 0.0:1.0:sum(map(maximum, (t1, t2, t3)))
+for i=1:size(corrm)[1], j=1:size(corrm)[2]
+    ti = tag(i)
+    tj = tag(j)
+    energy!(s, ti, ev2angphz(energies[i])-frame)
+    dipole!(s, "G", ti, tdm[i])
+    dipole!(s, ti, tj, 0)
+    lineshape!(s, "G", ti, zero)
+    lineshape!(s, ti, "G", zero)
+    #if i == j
+        lut = LineshapeLUT(t->(g_inhomo(t, ev2angphz(corrm[i,j]))
+                               +g_homo(t, ev2angphz(0.005))
+                               +g_inhomo(t, ev2angphz(sigma))), lut_grid)
+    #else
+    # we assume perfect correlation along the inhomogeneous coordinate.
+    #   lut = LineshapeLUT(t->(g_inhomo(t, ev2angphz(corrm[i,j]))
+    #                        +g_inhomo(t, ev2angphz(sigma))), lut_grid)
+    #end
+    lineshape!(s, ti, tj, lut)
+end
+info("    Took $(now()-t0_sys)")
+info("Number of order 1 Hilbert Paths: $(length(collect(hilbert_paths(s, 1))))")
+info("Number of order 3 Hilbert Paths: $(length(collect(hilbert_paths(s, 3))))")
+
+grd_lin = TimeGrid(t1)
+
+info("Computing linear responses separately.")
+t0_calc = now()
+out_root = "$(splitext(basename(states_fn))[1])_$(splitext(basename(cf_fn))[1])"
+
+# compute them separately
+totlin = zeros(Complex128, size(grd_lin))
+for p in hilbert_paths(s, 1)
+    rlin = linear(grd_lin, s, p)
+    totlin += rlin
+    out_lin = out_root*"_lin_"*join(p.p[2:end-1], "-")*".txt"
+    info("Saving linear to: $out_lin")
+    writedlm(out_lin, [grid(grd_lin)[1] real(rlin) imag(rlin)])
+end
+writedlm(out_root*"_lin_tot.txt", [grid(grd_lin)[1] real(totlin) imag(totlin)])
+info("    Took $(now()-t0_calc)")
+info("Preparing calculation of third order response")
+grd_trd = TimeGrid(t1, t2, t3)
+info("Grid size: $(size(grd_trd))")
+
+rr_fn = out_root*"_rr.bin"
+rn_fn = out_root*"_rn.bin"
+info("Saving rephasing to: $rr_fn")
+info("Saving nonrephasing to: $rr_fn")
+rr_out = open(rr_fn, "w+")
+rn_out = open(rn_fn, "w+")
+rr = Mmap.mmap(rr_out, Array{Complex128, 3}, size(grd_trd))
+rn = Mmap.mmap(rn_out, Array{Complex128, 3}, size(grd_trd))
+pm = Progress(4*length(collect(hilbert_paths(s, 3))))
+info("Computing third order response")
+t0_third = now()
+rr .+= R2(grd_trd, s, pm)
+rr .+= R3(grd_trd, s, pm)
+rn .+= R1(grd_trd, s, pm)
+rn .+= R4(grd_trd, s, pm)
+info("    Took $(now()-t0_third)")
+close(rr_out)
+close(rn_out)
+info("Total runtime: $(now()-t0)")
+end # run
+
+run(ARGS)
diff --git a/examples/aimd_static/xcorr_0_10.txt b/examples/aimd_static/xcorr_0_10.txt
@@ -0,0 +1,10 @@
+4.087452459290188272e-02 3.630889293716435540e-02 3.286830737009410752e-02 3.382167758160792292e-02 3.422752399138342949e-02 2.809560537724520649e-02 3.300191972712920863e-02 2.501988921409759384e-02 2.279167471652290139e-02 2.609973013632716290e-02
+3.630889293716435540e-02 3.548786907339005620e-02 2.990804695203031285e-02 2.830750100865434687e-02 2.867135074440373826e-02 2.642927010813047456e-02 3.002235141376016422e-02 2.316426700524965715e-02 2.010527027931844832e-02 2.247395102925969779e-02
+3.286830737009410752e-02 2.990804695203031285e-02 3.173017607196557816e-02 2.870322593045149945e-02 2.810510740025148890e-02 2.300875665305159717e-02 2.756556054995471655e-02 2.250985127726257834e-02 2.121887819665617056e-02 2.365889249231832928e-02
+3.382167758160792292e-02 2.830750100865434687e-02 2.870322593045149945e-02 3.285563114019410258e-02 3.121266978998517641e-02 2.308462075196793251e-02 2.788389502491780969e-02 2.203627535989740802e-02 2.188458636072978733e-02 2.476908182934962205e-02
+3.422752399138342949e-02 2.867135074440373826e-02 2.810510740025148890e-02 3.121266978998517641e-02 3.323494132043854771e-02 2.463587600880751821e-02 2.852101694886700206e-02 2.328042367987211841e-02 2.308167736070790982e-02 2.576519414747168585e-02
+2.809560537724520649e-02 2.642927010813047456e-02 2.300875665305159717e-02 2.308462075196793251e-02 2.463587600880751821e-02 2.356163794577150958e-02 2.385205758254578534e-02 1.935496707159090057e-02 1.834417230938027721e-02 1.979222504115265430e-02
+3.300191972712920863e-02 3.002235141376016422e-02 2.756556054995471655e-02 2.788389502491780969e-02 2.852101694886700206e-02 2.385205758254578534e-02 2.873707343835913558e-02 2.212929276798631176e-02 2.051240658578280807e-02 2.301343215528761643e-02
+2.501988921409759384e-02 2.316426700524965715e-02 2.250985127726257834e-02 2.203627535989740802e-02 2.328042367987211841e-02 1.935496707159090057e-02 2.212929276798631176e-02 2.023299591951278203e-02 2.040160044953351695e-02 2.042871080739329381e-02
+2.279167471652290139e-02 2.010527027931844832e-02 2.121887819665617056e-02 2.188458636072978733e-02 2.308167736070790982e-02 1.834417230938027721e-02 2.051240658578280807e-02 2.040160044953351695e-02 2.480529741786694289e-02 2.197884204119718646e-02
+2.609973013632716290e-02 2.247395102925969779e-02 2.365889249231832928e-02 2.476908182934962205e-02 2.576519414747168585e-02 1.979222504115265430e-02 2.301343215528761643e-02 2.042871080739329381e-02 2.197884204119718646e-02 2.290775574235751066e-02
