Skip to content

Commit

Permalink
exports now include zero step and docs have more function comments
Browse files Browse the repository at this point in the history
  • Loading branch information
tjira committed Jan 3, 2025
1 parent 09eee3a commit a4deb28
Show file tree
Hide file tree
Showing 17 changed files with 261 additions and 103 deletions.
149 changes: 77 additions & 72 deletions src/classicaldynamics.zig

Large diffs are not rendered by default.

10 changes: 10 additions & 0 deletions src/configurationinteraction.zig
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
//! Module for Configuration Interaction (CI) calculations.

const std = @import("std");

const hfm = @import("hartreefock.zig");
Expand All @@ -12,6 +14,7 @@ const Vector = @import("vector.zig").Vector;
const asfloat = @import("helper.zig").asfloat;
const c = @import("helper.zig").c ;

/// The CI options.
pub fn ConfigurationInteractionOptions(comptime T: type) type {
return struct {
excitation: ?[]const u32 = null,
Expand All @@ -20,16 +23,19 @@ pub fn ConfigurationInteractionOptions(comptime T: type) type {
};
}

/// The CI output.
pub fn ConfigurationInteractionOutput(comptime T: type) type {
return struct {
E: T,

/// Free the memory allocated for the CI output.
pub fn deinit(self: ConfigurationInteractionOutput(T)) void {
_ = self;
}
};
}

/// The main function to run the CI calculations.
pub fn run(comptime T: type, opt: ConfigurationInteractionOptions(T), print: bool, allocator: std.mem.Allocator) !ConfigurationInteractionOutput(T) {
const hf = try hfm.run(T, opt.hartree_fock, print, allocator);

Expand Down Expand Up @@ -90,6 +96,7 @@ pub fn run(comptime T: type, opt: ConfigurationInteractionOptions(T), print: boo
};
}

/// Aligns the vector C to the vector B. The result is stored in the vector A and the sign of the permutation is returned.
fn alignDeterminant(A: *Vector(usize), B: Vector(usize), C: Vector(usize)) !i32 {
@memcpy(A.data, C.data); var k: i32 = 0; var sign: i32 = 1;

Expand All @@ -102,6 +109,7 @@ fn alignDeterminant(A: *Vector(usize), B: Vector(usize), C: Vector(usize)) !i32
return sign;
}

/// Generates the determinants for the CI calculations.
fn generateDeterminants(D: *Matrix(usize), nbf: usize) void {
for (0..D.cols) |i| D.ptr(0, i).* = i;

Expand All @@ -117,6 +125,7 @@ fn generateDeterminants(D: *Matrix(usize), nbf: usize) void {
}
}

/// Slater-Condon rules for the CI calculations.
fn slater(comptime T: type, A: Vector(usize), so: []const usize, H_MS: Matrix(T), J_MS_A: Tensor(T)) !T {
var hij: T = 0;

Expand Down Expand Up @@ -147,6 +156,7 @@ fn slater(comptime T: type, A: Vector(usize), so: []const usize, H_MS: Matrix(T)
return hij;
}

/// Function to perform all integrals transformations used in the CI calculations.
fn transform(comptime T: type, H_MS: *Matrix(T), J_MS_A: *Tensor(T), T_AO: Matrix(T), V_AO: Matrix(T), J_AO: Tensor(T), C_MO: Matrix(T), allocator: std.mem.Allocator) !void {
var H_AO = try Matrix(T).init(T_AO.rows, T_AO.cols, allocator); defer H_AO.deinit(); mat.add(T, &H_AO, T_AO, V_AO);

Expand Down
4 changes: 4 additions & 0 deletions src/constant.zig
Original file line number Diff line number Diff line change
@@ -1,7 +1,11 @@
//! Constants used in the project.

const std = @import("std");

/// Angstrom to atomic units factor.
pub const A2AU = 1.8897261254578281;

/// Atomic symbol to atomic number map.
pub const SM2AN = std.StaticStringMap(u32).initComptime(.{
.{ "H", 1}, .{"He", 2},
.{"Li", 3}, .{"Be", 4}, .{ "B", 5}, .{ "C", 6}, .{ "N", 7}, .{ "O", 8}, .{ "F", 9}, .{"Ne", 10},
Expand Down
4 changes: 4 additions & 0 deletions src/fouriertransform.zig
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
//! Fourier transform module.

const std = @import("std"); const Complex = std.math.Complex;

const vec = @import("vector.zig");
Expand All @@ -9,6 +11,7 @@ const asfloat = @import("helper.zig").asfloat;
const bitrev = @import("helper.zig").bitrev ;
const prod = @import("helper.zig").prod ;

/// Fast Fourier transform for a one-dimensional array. The factor argument is the value in the exponent of the Fourier transform. Factor -1 corresponds to the forward Fourier transform, while factor 1 corresponds to the inverse Fourier transform.
pub fn fft(comptime T: type, arr: StridedArray(Complex(T)), factor: i32) !void {
const n = arr.len; const logn: u6 = @intCast(std.math.log2(n));

Expand Down Expand Up @@ -50,6 +53,7 @@ pub fn fft(comptime T: type, arr: StridedArray(Complex(T)), factor: i32) !void {
};
}

/// Fast Fourier transform for an n-dimensional array. The factor argument is the value in the exponent of the Fourier transform. Factor -1 corresponds to the forward Fourier transform, while factor 1 corresponds to the inverse Fourier transform.
pub fn fftn(comptime T: type, arr: []Complex(T), shape: []const usize, factor: i32) !void {
const sprod = prod(usize, shape); var stride: usize = 1;

Expand Down
7 changes: 7 additions & 0 deletions src/hartreefock.zig
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
//! Module to perform the Hartree-Fock calculation.

const std = @import("std");

const A2AU = @import("constant.zig").A2AU ;
Expand All @@ -14,6 +16,7 @@ const Tensor = @import("tensor.zig").Tensor;
const asfloat = @import("helper.zig").asfloat;
const uncr = @import("helper.zig").uncr ;

/// The Hartree-Fock options.
pub fn HartreeFockOptions(comptime T: type) type {
return struct {
const Integral = struct {
Expand All @@ -31,16 +34,19 @@ pub fn HartreeFockOptions(comptime T: type) type {
};
}

/// The Hartree-Fock output.
pub fn HartreeFockOutput(comptime T: type) type {
return struct {
C_MO: Matrix(T), D_MO: Matrix(T), E_MO: Matrix(T), F_AO: Matrix(T), E: T, VNN: T, nbf: usize, nocc: usize,

/// Free the memory allocated for the Hartree-Fock output.
pub fn deinit(self: HartreeFockOutput(T)) void {
self.C_MO.deinit(); self.D_MO.deinit(); self.E_MO.deinit(); self.F_AO.deinit();
}
};
}

/// Run the Hartree-Fock calculation with the given options.
pub fn run(comptime T: type, opt: HartreeFockOptions(T), print: bool, allocator: std.mem.Allocator) !HartreeFockOutput(T) {
const S_AO = try mat.read(T, opt.integral.overlap, allocator); defer S_AO.deinit();

Expand Down Expand Up @@ -118,6 +124,7 @@ pub fn run(comptime T: type, opt: HartreeFockOptions(T), print: bool, allocator:
};
}

/// Parse the .xyz system from the given path.
pub fn parseSystem(comptime T: type, path: []const u8, print: bool, allocator: std.mem.Allocator) !struct {natom: u32, nocc: u32, VNN: T} {
const file = try std.fs.cwd().openFile(path, .{}); defer file.close(); var buffer: [64]u8 = undefined;

Expand Down
7 changes: 7 additions & 0 deletions src/helper.zig
Original file line number Diff line number Diff line change
@@ -1,21 +1,28 @@
//! Module with some helper functions.

const std = @import("std");

/// Casts an integer to a float.
pub fn asfloat(comptime T: type, a: anytype) T {
return @as(T, @floatFromInt(a));
}

/// Reverse the bits of a number.
pub fn bitrev(value: anytype, count: u6) @TypeOf(value) {
var result: @TypeOf(value) = 0; var i: u6 = 0; while (i < count) : (i += 1) {result |= ((value >> @intCast(i)) & 1) << @intCast(count - 1 - i);} return result;
}

/// Calculate the combination number.
pub fn c(n: anytype, k: @TypeOf(n)) @TypeOf(n) {
var nck: @TypeOf(n) = 1; for (k + 1..n + 1) |i| {nck *= i;} for (2..n - k + 1) |i| {nck /= i;} return nck;
}

/// Calculate the product of an array.
pub fn prod(comptime T: type, v: []const T) T {
var result: T = 1; for (v) |value| result *= value; return result;
}

/// Remove the carriage return from a string. Fucking windows.
pub fn uncr(string: []const u8) []const u8 {
return if (string[string.len - 1] == 13) string[0..string.len - 1] else string;
}
5 changes: 1 addition & 4 deletions src/main.zig
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ pub const Matrix = @import("matrix.zig").Matrix;
pub const Tensor = @import("tensor.zig").Tensor;
pub const Vector = @import("vector.zig").Vector;

/// Parse the input JSON file and run the corresponding target.
pub fn parse(filebuf: []const u8) !void {
const inputjs = try std.json.parseFromSlice(std.json.Value, allocator, filebuf, .{}); defer inputjs.deinit();

Expand Down Expand Up @@ -91,7 +92,3 @@ pub fn main() !void {

std.debug.print("\nTOTAL EXECUTION TIME: {}\n", .{std.fmt.fmtDuration(timer.read())});
}

test {
std.testing.refAllDecls(@This());
}
28 changes: 28 additions & 0 deletions src/matrix.zig
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
//! This module provides a matrix class and functions to manipulate matrices.

const std = @import("std"); const Complex = std.math.Complex;

const StridedArray = @import("stridedarray.zig").StridedArray;
Expand All @@ -20,17 +22,20 @@ pub fn Matrix(comptime T: type) type {
.allocator = allocator
};
}
/// Free the memory allocated for the matrix.
pub fn deinit(self: Matrix(T)) void {
self.allocator.free(self.data);
}

/// Clone the matrix. The function returns an error if the allocation of the new matrix fails.
pub fn clone(self: Matrix(T)) !Matrix(T) {
const other = try Matrix(T).init(self.rows, self.cols, self.allocator);

@memcpy(other.data, self.data);

return other;
}
/// Convert the matrix to a complex matrix. The function returns an error if the allocation of the new matrix fails.
pub fn complex(self: Matrix(T)) !Matrix(Complex(T)) {
var other = try Matrix(Complex(T)).init(self.rows, self.cols, self.allocator);

Expand All @@ -41,12 +46,15 @@ pub fn Matrix(comptime T: type) type {
return other;
}

/// Access the element at the given row and column as a value.
pub fn at(self: Matrix(T), i: usize, j: usize) T {
return self.ptr(i, j).*;
}
/// Access the element at the given row and column as a mutable reference.
pub fn ptr(self: Matrix(T), i: usize, j: usize) *T {
return &self.data[i * self.cols + j];
}
/// Returns a reference to the row at the given index. The row is returned as a new matrix. No memory is allocated.
pub fn row(self: Matrix(T), i: usize) Matrix(T) {
return Matrix(T){
.data = self.data[i * self.cols..(i + 1) * self.cols],
Expand All @@ -55,9 +63,11 @@ pub fn Matrix(comptime T: type) type {
.allocator = self.allocator
};
}
/// Returns the matrix in the form of a strided array. No memory is allocated.
pub fn sa(self: Matrix(T)) StridedArray(T) {
return StridedArray(T){.data = self.data, .len = self.rows * self.cols, .stride = 1, .zero = 0};
}
/// Returns the matrix in the form of a vector. No memory is allocated.
pub fn vector(self: Matrix(T)) Vector(T) {
return Vector(T){
.data = self.data[0..],
Expand All @@ -66,29 +76,34 @@ pub fn Matrix(comptime T: type) type {
};
}

/// Fill the matrix with a given value.
pub fn fill(self: Matrix(T), value: T) void {
for (0..self.data.len) |i| self.data[i] = value;
}
/// Fill the diagonal of the matrix with ones.
pub fn identity(self: Matrix(T)) void {
self.fill(0);

for (0..self.rows) |i| {
self.ptr(i, i).* = 1;
}
}
/// Fill the data of the matrix with a linearly spaced vector from start to end. Both start and end are included.
pub fn linspace(self: Matrix(T), start: T, end: T) void {
for (0..self.data.len) |i| {
self.data[i] = start + asfloat(T, i) * (end - start) / asfloat(T, self.rows * self.cols - 1);
}
}

/// Print the matrix to the given device.
pub fn print(self: Matrix(T), device: anytype) !void {
try device.print("{d} {d}\n", .{self.rows, self.cols});

for (self.data, 1..) |e, i| {
try device.print("{d:20.14}{s}", .{e, if(i % self.cols == 0) "\n" else " "});
}
}
/// Write the matrix to a file.
pub fn write(self: Matrix(T), path: []const u8) !void {
const file = try std.fs.cwd().createFile(path, .{}); defer file.close();

Expand All @@ -97,6 +112,7 @@ pub fn Matrix(comptime T: type) type {
};
}

/// Add two matrices element-wise. The output matrix is stored in the matrix C.
pub fn add(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
if (@typeInfo(T) != .Struct) for (0..C.data.len) |i| {
C.data[i] = A.data[i] + B.data[i];
Expand All @@ -107,6 +123,7 @@ pub fn add(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
};
}

/// Add a scalar to a matrix element-wise. The output matrix is stored in the matrix C.
pub fn adds(comptime T: type, C: *Matrix(T), A: Matrix(T), s: T) void {
if (@typeInfo(T) != .Struct) for (0..C.data.len) |i| {
C.data[i] = A.data[i] + s;
Expand All @@ -117,6 +134,7 @@ pub fn adds(comptime T: type, C: *Matrix(T), A: Matrix(T), s: T) void {
};
}

/// Find the eigenvalues and eigenvectors of a real symmetric matrix A. The eigenvalues are stored in the diagonal of the matrix J, and the eigenvectors are stored in the columns of the matrix C. The matrices T1 and T2 are temporary matrices used in the computation.
pub fn eigh(comptime T: type, J: *Matrix(T), C: *Matrix(T), A: Matrix(T), T1: *Matrix(T), T2: *Matrix(T)) void {
var maxi: usize = 0; var maxj: usize = 1; var maxv: T = 0; var phi: T = undefined; @memcpy(J.data, A.data); C.identity();

Expand Down Expand Up @@ -150,6 +168,7 @@ pub fn eigh(comptime T: type, J: *Matrix(T), C: *Matrix(T), A: Matrix(T), T1: *M
};
}

/// Check if two matrices are equal within a given tolerance. The function returns true if the matrices are equal and false otherwise.
pub fn eq(comptime T: type, A: Matrix(T), B: Matrix(T), epsilon: T) bool {
if (A.rows != B.rows or A.cols != B.cols) return false;

Expand All @@ -164,6 +183,7 @@ pub fn eq(comptime T: type, A: Matrix(T), B: Matrix(T), epsilon: T) bool {
return true;
}

/// Horizontally concatenate two matrices. The output matrix is stored in the matrix C.
pub fn hjoin(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
for (0..A.rows) |i| {

Expand All @@ -177,6 +197,7 @@ pub fn hjoin(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
}
}

/// Matrix multiplication of the adjoint of matrix A and matrix B. The output matrix is stored in the matrix C.
pub fn mam(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
if (@typeInfo(T) == .Struct) C.fill(T.init(0, 0)) else C.fill(0);

Expand All @@ -189,6 +210,7 @@ pub fn mam(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
};
}

/// Matrix multiplication of matrix A and matrix B. The output matrix is stored in the matrix C.
pub fn mm(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
if (@typeInfo(T) == .Struct) C.fill(T.init(0, 0)) else C.fill(0);

Expand All @@ -201,6 +223,7 @@ pub fn mm(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
};
}

/// Matrix multiplication of matrix A and the adjoint of matrix B. The output matrix is stored in the matrix C.
pub fn mma(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
if (@typeInfo(T) == .Struct) C.fill(T.init(0, 0)) else C.fill(0);

Expand All @@ -213,6 +236,7 @@ pub fn mma(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
};
}

/// Multiply two matrices element-wise. The output matrix is stored in the matrix C.
pub fn mul(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
if (@typeInfo(T) != .Struct) for (0..C.data.len) |i| {
C.data[i] = A.data[i] * B.data[i];
Expand All @@ -223,6 +247,7 @@ pub fn mul(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
};
}

/// Multiply a matrix by a scalar element-wise. The output matrix is stored in the matrix C.
pub fn muls(comptime T: type, C: *Matrix(T), A: Matrix(T), s: T) void {
if (@typeInfo(T) != .Struct) for (0..C.data.len) |i| {
C.data[i] = A.data[i] * s;
Expand All @@ -233,6 +258,7 @@ pub fn muls(comptime T: type, C: *Matrix(T), A: Matrix(T), s: T) void {
};
}

/// Read a matrix from a file. The function returns an error if the file cannot be opened or if the matrix cannot be read.
pub fn read(comptime T: type, path: []const u8, allocator: std.mem.Allocator) !Matrix(T) {
const file = try std.fs.cwd().openFile(path, .{}); defer file.close(); var buffer: [16]u8 = undefined;

Expand All @@ -256,6 +282,7 @@ pub fn read(comptime T: type, path: []const u8, allocator: std.mem.Allocator) !M
return mat;
}

/// Subtract two matrices element-wise. The output matrix is stored in the matrix C.
pub fn sub(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
if (@typeInfo(T) != .Struct) for (0..C.data.len) |i| {
C.data[i] = A.data[i] - B.data[i];
Expand All @@ -266,6 +293,7 @@ pub fn sub(comptime T: type, C: *Matrix(T), A: Matrix(T), B: Matrix(T)) void {
};
}

/// Transpose the matrix A. The output matrix is stored in the matrix C.
pub fn transpose(comptime T: type, C: *Matrix(T), A: Matrix(T)) void {
for (0..A.rows) |i| for (0..A.cols) |j| {
C.ptr(j, i).* = A.at(i, j);
Expand Down
Loading

0 comments on commit a4deb28

Please sign in to comment.