
// Robert Brunner
// 2/25/2000
//
// A small program to solve the dynamics of a system of atom-like
// particles with random charges and positions in a box.
// The particles interact through electrostatic interaction,
// with a cutoff.
// Problems:
//   - The particles may be randomly created very close to each other
//     causing huge energies and errors
//   - Since a cell-based approach is used, but atoms are not moved
//     to their correct cells as the simulation procedes, eventually
//     interacting pairs will be ignored.  To fix this, the atoms should
//     be frequently put back into their correct cells as they move in
//     space. (For this assignment, we will ignore this).
//   - In the force computation, and in integration, efficiency is
//     sometimes sacrificed for clarity and ease of debugging.

#define NEW_UPDATE_CELL

struct Atom {
  double m;     // Mass in AMU (= 1.66E-27 kg)
  double q;     // Charge in units of elemental charge ( = 1.602E-19C)
  double x, y, z;     // Position in Angstroms ( = 1E-10 m)
  double vx, vy, vz;    // Velocity in angstroms/fs ( = 1E5m/s)
  double vhx, vhy, vhz;  // Half-step velocity, used for Verlet integration
  double fx, fy, fz;    // Force in newtons
};

struct Cell {
  int n_atoms;
  int my_x, my_y, my_z;  // index of the cell
  double min_x, min_y, min_z;
  double max_x, max_y, max_z;   // probably unnecessary?
  Atom* atoms;   // array to be allocated ;
};

struct SimParams {
  int steps;
  int n_atoms;
  double min_x, min_y, min_z;
  double max_x, max_y, max_z;  // Angstroms
  int cell_dim_x, cell_dim_y, cell_dim_z;
  double cutoff;  // Angstroms
  double margin;  // Angstroms
  double step_sz; // femtoseconds (=1E-15s)
};

// A few global variables 
extern SimParams params;
extern double pot_energy;
extern double kin_energy;
extern bool first_step;

// Functions in littleMD.C
void init(void);
void done(void);
void calc_forces(void);
void update_coordinates(void);
 
// Functions in utilities.C
void init_cell(Cell** cell, int xi, int yi, int zi);
double cell_self(Cell* this_cell);
double cell_neighbor(Cell* cell1, Cell* cell2);
double calc_force(Atom *atom1 , Atom* atom2);
#ifdef NEW_UPDATE_CELL
double update_cell(Cell* this_cell, int first_step);
#else
double update_cell(Cell* this_cell);
#endif
