Changes

milestones/01/main.cpp

@@ -7,7 +7,6 @@
 #include <mpi.h>
 #endif
 
-
 int main(int argc, char *argv[]) {
     int rank = 0, size = 1;
 

milestones/01/meson.build

@@ -4,4 +4,4 @@ executable(
 	include_directories : [lib_incdirs],
 	link_with : [lib],
 	dependencies : [eigen, mpi]
-)
+)

milestones/02/main.cpp

@@ -0,0 +1,32 @@
+//
+// Created by 13396 on 2024/5/25.
+//
+#include <Eigen/Dense>
+#include <iostream>
+#include <verlet.h>
+
+int main() {
+
+  int nb_steps = 10;
+  double x, y, z, vx, vy, vz = 0;
+  double fx, fy ,fz = 1;
+  double timestep = 0.5;
+  double mass = 0.2;
+
+  for (int i = 0; i < nb_steps; ++i) {
+    double vx1 = vx;
+    double vy1 = vy;
+    double vz1 = vz;
+    std::cout << "Step: " << i << std::endl;
+
+    verlet_step1(x, y, z, vx, vy, vz, fx, fy, fz, timestep, mass);
+
+    fx = mass * (vx - vx1)/ timestep;
+    fy = mass * (vy - vy1)/ timestep;
+    fz = mass * (vz - vz1)/ timestep;
+
+    verlet_step2(vx, vy, vz, fx, fy, fz, timestep, mass);
+      std::cout << "x" << x << std::endl;
+  }
+  return 0;
+}

milestones/02/meson.build

@@ -0,0 +1,7 @@
+executable(
+    'milestone02',
+    'main.cpp',
+    include_directories : [lib_incdirs],
+    link_with : [lib],
+    dependencies : [eigen, mpi]
+)

milestones/04/main.cpp

@@ -0,0 +1,73 @@
+//
+// Created by 13396 on 2024/5/25.
+//
+
+#include <iostream>
+#include <lj_direction_summation.cpp>
+
+#include <string>
+#include <xyz.h>
+
+double compute_kinetic_energy(const Eigen::Array3Xd &velocities, double mass) {
+    double kinetic_energy = 0.0;
+    size_t num_atoms = velocities.cols();
+    for (int i = 0; i < num_atoms; ++i) {
+        kinetic_energy += 0.5 * mass * velocities.col(i).matrix().squaredNorm();
+    }
+    return kinetic_energy;
+}
+
+void velocity_verlet(Atoms &atoms, double dt, double epsilon, double sigma, double mass) {
+    size_t num_atoms = atoms.nb_atoms();
+    // Velocity Verlet Integration
+    for (int i = 0; i < num_atoms; ++i) {
+        // Update positions
+        atoms.positions.col(i) += atoms.velocities.col(i) * dt + 0.5 * atoms.forces.col(i) * dt * dt / mass;
+    }
+
+    for (int i = 0; i < num_atoms; ++i) {
+        // Update velocities
+        atoms.velocities.col(i) += 0.5 * dt * (atoms.forces.col(i) + atoms.forces.col(i)) / mass;
+    }
+}
+
+
+int main() {
+try{
+    std::string file_path = "./lj54.xyz";
+    auto [names, positions, velocities]{read_xyz_with_velocities(file_path)};
+    Atoms atoms(positions, velocities);
+
+    double epsilon = 1.0;
+    double sigma = 1.0;
+    double mass = 1.0;
+    double dt = 0.001;
+    double total_time = 10.0;
+    int num_steps = static_cast<int>(total_time / dt);
+    int output_interval = 100;
+
+    std::ofstream traj("./traj.xyz",std::ios::out);
+
+    for (int step = 0; step < num_steps; ++step) {
+        if (step % output_interval == 0) {
+            write_xyz(traj, atoms);
+        }
+
+        double kinetic_energy = compute_kinetic_energy(atoms.velocities, mass);
+        double potential_energy = lj_direct_summation(atoms, epsilon, sigma);
+        double total_energy = kinetic_energy + potential_energy;
+
+        std::cout << "Step: " << step
+                  << " Total Energy: " << total_energy
+                  << " Kinetic Energy: " << kinetic_energy
+                  << " Potential Energy: " << potential_energy << std::endl;
+
+        velocity_verlet(atoms, dt, epsilon, sigma, mass);
+    }
+    traj.close();
+} catch (const std::runtime_error &e) {
+        std::cerr << "Error: " << e.what() << std::endl;
+    }
+
+    return 0;
+}

milestones/04/meson.build

@@ -0,0 +1,22 @@
+executable(
+	'milestone04',
+	'main.cpp',
+	include_directories : [lib_incdirs],
+	link_with : [lib],
+	dependencies : [eigen, mpi]
+
+)
+
+#fs = import('fs')
+#fs.copyfile('lj54.xyz')
+
+# Define the source and destination paths
+source_file = 'lj54.xyz'
+destination_file = 'lj54.xyz'  # This is the name it will have in the build directory
+
+# Use configure_file to copy the file
+configure_file(
+	input : source_file,
+	output : destination_file,
+	copy: true
+)

milestones/05/main.cpp

@@ -0,0 +1,155 @@
+//
+// Created by 13396 on 2024/5/25.
+//
+#include <iostream>
+#include <lj_direction_summation.cpp>
+#include <string>
+#include <xyz.h>
+#include <chrono>
+
+double compute_kinetic_energy(const Eigen::Array3Xd &velocities, double mass) {
+    double kinetic_energy = 0.0;
+    size_t num_atoms = velocities.cols();
+    for (int i = 0; i < num_atoms; ++i) {
+        kinetic_energy += 0.5 * mass * velocities.col(i).matrix().squaredNorm();
+    }
+    return kinetic_energy;
+}
+
+void velocity_verlet(Atoms &atoms, double timestep, double epsilon, double sigma, double mass) {
+    size_t num_atoms = atoms.nb_atoms();
+    // Compute initial forces and potential energy
+    double potential_energy = lj_direct_summation(atoms, epsilon, sigma);
+
+    // Update positions
+    atoms.velocities += 0.5 * timestep * atoms.forces / mass; // Assuming mass = 1.0
+    atoms.positions += timestep * atoms.velocities;
+
+    // Compute forces and potential energy
+    potential_energy = lj_direct_summation(atoms, epsilon, sigma);
+
+    // Update velocities
+    atoms.velocities += 0.5 * timestep * atoms.forces / mass;
+}
+
+
+void apply_berendsen_thermostat(Atoms &atoms, double mass, double target_temp, double tau_T, double timestep) {
+    double current_temp = (2.0 / 3.0) * (compute_kinetic_energy(atoms.velocities, mass) / atoms.nb_atoms()); //(atoms.velocities.square().sum()) / (3.0 * atoms.nb_atoms());
+    double lambda = sqrt(1.0 + (timestep / tau_T) * (target_temp / current_temp - 1.0));
+    atoms.velocities *= lambda;
+}
+
+// Function to run the simulation
+double run_simulation(Atoms &atoms, double epsilon, double sigma, double mass,
+                    double target_temp, double tau_T, double timestep, double total_time,
+                    const std::string &output_file) {
+    int steps = static_cast<int>(total_time / timestep);
+
+    // File to record data
+    std::ofstream data(output_file);
+    data << "step,total_energy,kinetic_energy,potential_energy,temperature\n";
+
+    auto start = std::chrono::high_resolution_clock::now();
+
+    // Main MD loop
+    for (int step = 0; step <= steps; ++step) {
+        // Compute potential energy and forces
+        double potential_energy = lj_direct_summation(atoms, epsilon, sigma);
+
+        // Compute kinetic energy and temperature
+        double kinetic_energy = compute_kinetic_energy(atoms.velocities,mass);
+        double total_energy = potential_energy + kinetic_energy;
+        double temperature = (2.0 / 3.0) * (kinetic_energy / atoms.nb_atoms());
+
+        // Record data
+        data << step << "," << total_energy << "," << kinetic_energy << ","
+             << potential_energy << "," << temperature << "\n";
+
+        // Velocity Verlet integration
+        velocity_verlet(atoms, timestep, epsilon, sigma, mass);
+
+        // Apply Berendsen thermostat
+        apply_berendsen_thermostat(atoms, mass, target_temp, tau_T, timestep);
+        // std::cout << "velocity: " << atoms.velocities.col(1) << std::endl;
+    }
+
+    //Caculate the simulation time
+    auto end = std::chrono::high_resolution_clock::now();
+    std::chrono::duration<double> duration = end - start;
+    std::cout << "Number of atoms: " << atoms.nb_atoms() << ", Simulation time: " << duration.count() << " seconds" << std::endl;
+    double elapsed_time = duration.count();
+
+   //Close the write of data
+    data.close();
+
+    return elapsed_time;
+}
+
+void create_cubic_lattice(Atoms &atoms, int num_atoms_per_side, double lattice_constant) {
+    int num_atoms = num_atoms_per_side * num_atoms_per_side * num_atoms_per_side;
+    Positions_t positions(3, num_atoms);
+    Velocities_t velocities(3, num_atoms);
+    positions.setZero();
+    velocities.setZero();
+
+    int index = 0;
+    for (int x = 0; x < num_atoms_per_side; ++x) {
+        for (int y = 0; y < num_atoms_per_side; ++y) {
+            for (int z = 0; z < num_atoms_per_side; ++z) {
+                positions(0, index) = x * lattice_constant;
+                positions(1, index) = y * lattice_constant;
+                positions(2, index) = z * lattice_constant;
+                ++index;
+            }
+        }
+    }
+    atoms = Atoms(positions, velocities);
+}
+
+int main() {
+try{
+    //Parameters of atoms
+    Positions_t positions;
+    Velocities_t velocities;
+    Atoms atoms(positions, velocities);
+    double epsilon = 1.0;
+    double sigma = 1.0;
+    double mass = 1.0;
+
+    //Parameters of time
+    double timestep = 0.001;
+    double total_time = 2.5;
+    double target_temp = 1;  // Target temperature for the thermostat
+    double tau_T = 1.0;  // Thermostat coupling constant
+
+    //Marix of amount of atoms and simulation time
+    std::vector<int> atom_counts = {8, 27, 64, 125, 216, 343, 512, 729, 1000};
+    std::vector<double> times;
+
+    //Simulation with different numbers of atoms
+    for (int num_atoms_per_side = 2; num_atoms_per_side <= 10; ++num_atoms_per_side) {
+        double lattice_constant = 1.0;
+        create_cubic_lattice(atoms, num_atoms_per_side, lattice_constant);
+
+        std::string output_file = "simulation_data_" + std::to_string(atoms.nb_atoms()) + ".csv";
+
+        double elapsed_time = run_simulation(atoms, epsilon, sigma, mass,
+                            target_temp, tau_T, timestep,  total_time,
+                            output_file);
+
+        times.push_back(elapsed_time);
+    }
+
+    // Save the results to a file
+    std::ofstream result_file("simulation_times.csv");
+    result_file << "num_atoms,elapsed_time\n";
+    for (int i = 0; i < atom_counts.size(); ++i) {
+        result_file << atom_counts[i] << "," << times[i] << "\n";
+    }
+    result_file.close();
+} catch (const std::runtime_error &e) {
+        std::cerr << "Error: " << e.what() << std::endl;
+    }
+
+    return 0;
+}

milestones/05/meson.build

@@ -0,0 +1,19 @@
+executable(
+	'milestone05',
+	'main.cpp',
+	include_directories : [lib_incdirs],
+	link_with : [lib],
+	dependencies : [eigen, mpi]
+
+)
+
+# Define the source and destination paths
+source_file = 'lj54.xyz'
+destination_file = 'lj54.xyz'  # This is the name it will have in the build directory
+
+# Use configure_file to copy the file
+configure_file(
+	input : source_file,
+	output : destination_file,
+	copy: true
+)