change defaults for classical and quantum dynamics

tjira · Jan 15, 2025 · 6505c2c · 6505c2c
1 parent b2a86f2
commit 6505c2c
Show file tree

Hide file tree

Showing 7 changed files with 61 additions and 406 deletions.
diff --git a/.gitignore b/.gitignore
@@ -12,9 +12,9 @@ docs/**/*.pdf*
 
 education/**/*.mat*
 
-research/fssh_vs_lzsh/*.jso*
-research/fssh_vs_lzsh/*.mat*
-research/fssh_vs_lzsh/*.png*
+research/fssh_vs_lzsh/**/*.jso*
+research/fssh_vs_lzsh/**/*.mat*
+research/fssh_vs_lzsh/**/*.png*
 
 research/uracil_lvc/*.mat*
 research/uracil_lvc/*.png*

diff --git a/research/fssh_vs_lzsh/mash_check/exact.json b/research/fssh_vs_lzsh/mash_check/exact.json
@@ -1,7 +1,7 @@
 {
     "quantum_dynamics" : {
         "adiabatic" : true,
-        "imaginary" : false,
+        "mode" : [0, 1],
         "iterations" : 350,
         "time_step" : 10,
         "grid" : {

diff --git a/research/fssh_vs_lzsh/mash_check/mash.json b/research/fssh_vs_lzsh/mash_check/mash.json
@@ -31,7 +31,7 @@
         "potential" : "tully1D_3",
         "spin_mapping" : {
             "quantum_substep" : 10,
-            "sample_bloch_vector" : true
+            "sample_bloch_vector" : false
         },
         "time_derivative_coupling" : "numeric"
     }

diff --git a/src/classicaldynamics.zig b/src/classicaldynamics.zig
@@ -14,11 +14,11 @@ const asfloat = @import("helper.zig").asfloat;
 pub fn ClassicalDynamicsOptions(comptime T: type) type {
     return struct {
         const InitialConditions = struct {
-            position_mean: []const T = &[_]T{-10.0},
-            position_std:  []const T = &[_]T{  0.5},
-            momentum_mean: []const T = &[_]T{ 15.0},
-            momentum_std:  []const T = &[_]T{  1.0},
-            state: u32 = 1, mass: T = 2000
+            position_mean: []const T = &[_]T{1.0},
+            position_std:  []const T = &[_]T{0.5},
+            momentum_mean: []const T = &[_]T{0.0},
+            momentum_std:  []const T = &[_]T{1.0},
+            state: u32 = 0, mass: T = 1
         };
         const FewestSwitches = struct {
             quantum_substep: u32 = 10,
@@ -46,12 +46,12 @@ pub fn ClassicalDynamicsOptions(comptime T: type) type {
 
         adiabatic: bool = true,
         derivative_step: T = 0.001,
-        iterations: u32 = 3000,
-        potential: []const u8 = "tully1D_1",
+        iterations: u32 = 100,
+        potential: []const u8 = "harmonic1D_1",
         seed: u32 = 1,
         time_derivative_coupling: ?[]const u8 = "numeric",
-        time_step: T = 1,
-        trajectories: u32 = 100,
+        time_step: T = 0.1,
+        trajectories: u32 = 1000,
 
         fewest_switches: ?FewestSwitches = null, landau_zener: ?LandauZener = null, spin_mapping: ?SpinMapping = null,
 
@@ -64,15 +64,27 @@ pub fn ClassicalDynamicsOutput(comptime T: type) type {
     return struct {
         pop: Vector(T),
 
+        r: Vector(T),
+        p: Vector(T),
+
+        Ekin: T,
+        Epot: T,
+
         /// Initialize the classical dynamics output.
-        pub fn init(nstate: usize, allocator: std.mem.Allocator) !ClassicalDynamicsOutput(T) {
+        pub fn init(ndim: usize, nstate: usize, allocator: std.mem.Allocator) !ClassicalDynamicsOutput(T) {
             return ClassicalDynamicsOutput(T){
-                .pop = try Vector(T).init(nstate, allocator)
+                .pop = try Vector(T).init(nstate, allocator),
+
+                .r = try Vector(T).init(ndim, allocator),
+                .p = try Vector(T).init(ndim, allocator),
+
+                .Ekin = 0,
+                .Epot = 0
             };
         }
         /// Free the memory allocated for the classical dynamics output.
         pub fn deinit(self: ClassicalDynamicsOutput(T)) void {
-            self.pop.deinit();
+            self.pop.deinit(); self.r.deinit(); self.p.deinit();
         }
     };
 }
@@ -83,7 +95,7 @@ pub fn run(comptime T: type, opt: ClassicalDynamicsOptions(T), print: bool, allo
 
     const ndim = try mpt.dims(opt.potential); const nstate = try mpt.states(opt.potential);
 
-    var output = try ClassicalDynamicsOutput(T).init(nstate, allocator);
+    var output = try ClassicalDynamicsOutput(T).init(ndim, nstate, allocator);
 
     var prng_jump = std.Random.DefaultPrng.init(opt.seed); const rand_jump = prng_jump.random();
     var prng_traj = std.Random.DefaultPrng.init(opt.seed); const rand_traj = prng_traj.random();
@@ -219,11 +231,19 @@ pub fn run(comptime T: type, opt: ClassicalDynamicsOptions(T), print: bool, allo
                 if (opt.write.total_energy_mean             != null) etot.ptr(j, 1 + 0).* += Ekin + Epot                                                     ;
                 if (opt.write.position_mean                 != null) for (0..r.rows) |k| {position.ptr(j, 1 + k).* += r.at(k);}                              ;
                 if (opt.write.momentum_mean                 != null) for (0..v.rows) |k| {momentum.ptr(j, 1 + k).* += v.at(k) * opt.initial_conditions.mass;};
-                if (opt.write.time_derivative_coupling_mean != null) for (0..tdc.cols) |k| {tdc.ptr(j, 1 + k).* += TDC.data[k];}                             ;
+                if (opt.write.time_derivative_coupling_mean != null) for (0..TDC.rows * TDC.cols) |k| {tdc.ptr(j, 1 + k).* += TDC.data[k];}                  ;
                 if (opt.write.fssh_coefficient_mean         != null) for (0..C.rows) |k| {coef.ptr(j, 1 + k).* += C.at(k).magnitude() * C.at(k).magnitude();};
 
                 if (j == opt.iterations - 1) {
+
                     output.pop.ptr(s).* += 1;
+
+                    for (0..ndim) |k| {
+                        output.r.ptr(k).* += r.at(k);  output.p.ptr(k).* += v.at(k) * opt.initial_conditions.mass;
+                    }
+
+                    output.Ekin += Ekin;
+                    output.Epot += Epot;
                 }
 
                 if (print and (i == 0 or (i + 1) % opt.log_intervals.trajectory == 0) and (j % opt.log_intervals.iteration == 0)) {
@@ -243,6 +263,11 @@ pub fn run(comptime T: type, opt: ClassicalDynamicsOptions(T), print: bool, allo
     for (0..opt.iterations + 1) |i| {for (1..coef.cols    ) |j|     coef.ptr(i, j).* /= asfloat(T, opt.trajectories);}
 
     for (0..output.pop.rows) |i| output.pop.ptr(i).* /= asfloat(T, opt.trajectories);
+    for (0..output.r.rows  ) |i|   output.r.ptr(i).* /= asfloat(T, opt.trajectories);
+    for (0..output.p.rows  ) |i|   output.p.ptr(i).* /= asfloat(T, opt.trajectories);
+
+    output.Ekin /= asfloat(T, opt.trajectories);
+    output.Epot /= asfloat(T, opt.trajectories);
 
     for (0..nstate) |i| {
         if (print) {try std.io.getStdOut().writer().print("{s}FINAL POPULATION OF STATE {d:2}: {d:.6}\n", .{if (i == 0) "\n" else "", i, output.pop.at(i)});}

diff --git a/src/quantumdynamics.zig b/src/quantumdynamics.zig
@@ -17,10 +17,10 @@ const asfloat = @import("helper.zig").asfloat;
 pub fn QuantumDynamicsOptions(comptime T: type) type {
     return struct {
         const Grid = struct {
-            limits: []const T = &[_]f64{-16, 32}, points: u32 = 512
+            limits: []const T = &[_]f64{-8, 8}, points: u32 = 32
         };
         const InitialConditions = struct {
-            position: []const T = &[_]f64{-10}, momentum: []const T = &[_]f64{15}, gamma: T = 2, state: u32 = 1, mass: T = 2000
+            position: []const T = &[_]f64{1}, momentum: []const T = &[_]f64{0}, gamma: T = 2, state: u32 = 0, mass: T = 1
         };
         const LogIntervals = struct {
             iteration: u32 = 1
@@ -36,11 +36,11 @@ pub fn QuantumDynamicsOptions(comptime T: type) type {
             wavefunction: ?[]const u8 = null
         };
 
-        adiabatic: bool = true,
-        iterations: u32 = 300,
+        adiabatic: bool = false,
+        iterations: u32 = 100,
         mode: []const u32 = &[_]u32{0, 1},
-        potential: []const u8 = "tully1D_1",
-        time_step: T = 10,
+        potential: []const u8 = "harmonic1D_1",
+        time_step: T = 0.1,
 
         grid: Grid = .{}, initial_conditions: InitialConditions = .{}, log_intervals: LogIntervals = .{}, write: Write = .{}
     };
@@ -52,17 +52,20 @@ pub fn QuantumDynamicsOutput(comptime T: type) type {
         P: Matrix(T),
         r: Vector(T),
         p: Vector(T),
+
         Ekin: T,
         Epot: T,
 
         /// Initialize the quantum dynamics output struct.
         pub fn init(ndim: usize, nstate: usize, allocator: std.mem.Allocator) !QuantumDynamicsOutput(T) {
             return QuantumDynamicsOutput(T){
                 .P = try Matrix(T).init(nstate, nstate, allocator),
+
                 .r = try Vector(T).init(ndim, allocator),
                 .p = try Vector(T).init(ndim, allocator),
-                .Ekin = undefined,
-                .Epot = undefined
+
+                .Ekin = 0,
+                .Epot = 0
             };
         }
         /// Free the memory allocated for the quantum dynamics output struct.

diff --git a/test/classicaldynamics.zig b/test/classicaldynamics.zig
@@ -5,194 +5,10 @@ const vec = @import("acorn").vec;
 
 const Vector = @import("acorn").Vector;
 
-test "cdyn_fssh_doubleState1D_1" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.70; pop.ptr(1).* = 0.30;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "doubleState1D_1"; opt.fewest_switches = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_doubleState1D_1" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.67; pop.ptr(1).* = 0.33;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "doubleState1D_1"; opt.landau_zener = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_fssh_doubleState1D_2" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.96; pop.ptr(1).* = 0.04;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "doubleState1D_2"; opt.fewest_switches = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_doubleState1D_2" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.92; pop.ptr(1).* = 0.08;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "doubleState1D_2"; opt.landau_zener = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_fssh_tripleState1D_1" {
-    var pop = try Vector(f64).init(3, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.49; pop.ptr(1).* = 0.45; pop.ptr(2).* = 0.06;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tripleState1D_1"; opt.fewest_switches = .{}; opt.initial_conditions.state = 2;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_tripleState1D_1" {
-    var pop = try Vector(f64).init(3, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.41; pop.ptr(1).* = 0.59; pop.ptr(2).* = 0.00;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tripleState1D_1"; opt.landau_zener = .{}; opt.initial_conditions.state = 2;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_fssh_tripleState1D_2" {
-    var pop = try Vector(f64).init(3, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.96; pop.ptr(1).* = 0.04; pop.ptr(2).* = 0.00;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tripleState1D_2"; opt.fewest_switches = .{}; opt.initial_conditions.state = 2;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_tripleState1D_2" {
-    var pop = try Vector(f64).init(3, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.47; pop.ptr(1).* = 0.53; pop.ptr(2).* = 0.00;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tripleState1D_2"; opt.landau_zener = .{}; opt.initial_conditions.state = 2;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_fssh_tripleState1D_3" {
-    var pop = try Vector(f64).init(3, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.13; pop.ptr(1).* = 0.79; pop.ptr(2).* = 0.08;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tripleState1D_3"; opt.fewest_switches = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_tripleState1D_3" {
-    var pop = try Vector(f64).init(3, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.19; pop.ptr(1).* = 0.75; pop.ptr(2).* = 0.06;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tripleState1D_3"; opt.landau_zener = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_fssh_tully1D_1" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.40; pop.ptr(1).* = 0.60;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tully1D_1"; opt.fewest_switches = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_tully1D_1" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.51; pop.ptr(1).* = 0.49;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tully1D_1"; opt.landau_zener = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_fssh_tully1D_2" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.11; pop.ptr(1).* = 0.89;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tully1D_2"; opt.fewest_switches = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_tully1D_2" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.07; pop.ptr(1).* = 0.93;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tully1D_2"; opt.landau_zener = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_fssh_tully1D_3" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0.59; pop.ptr(1).* = 0.41;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tully1D_3"; opt.fewest_switches = .{}; opt.initial_conditions.state = 1;
-
-    const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
-
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
-}
-
-test "cdyn_lzsh_tully1D_3" {
-    var pop = try Vector(f64).init(2, std.testing.allocator); defer pop.deinit();
-
-    pop.ptr(0).* = 0; pop.ptr(1).* = 1;
-
-    var opt = cdn.ClassicalDynamicsOptions(f64){}; opt.potential = "tully1D_3"; opt.landau_zener = .{}; opt.initial_conditions.state = 1;
+test "cdyn_default" {
+    const opt = cdn.ClassicalDynamicsOptions(f64){};
 
     const output = try cdn.run(f64, opt, false, std.testing.allocator); defer output.deinit();
 
-    try std.testing.expect(vec.eq(f64, output.pop, pop, 1e-12));
+    try std.testing.expect(@abs(output.Ekin + output.Epot - 1.09553019324015) < 1e-8);
 }