/*
 * Hacking atoms
 * Units: mass in 10^-27 kg, distance in nanometers, time in picoseconds,
 * force in piconewtons, energy in milli-attojoules (10^-21 J)
 */

#include "atom.h"

void atom::
iterate_motion (double time_step)
{
  double m, a0, a1, a2;
  
  m = 0.5 * time_step / mass;
  
  a0 = m * f[0];
  a1 = m * f[1];
  a2 = m * f[2];
  
  v[0] += a0;
  v[1] += a1;
  v[2] += a2;
  
  x[0] += v[0] * time_step;
  x[1] += v[1] * time_step;
  x[2] += v[2] * time_step;
  
  v[0] += a0;
  v[1] += a1;
  v[2] += a2;
}

/* it might seem silly, but it helps keep main() a little simpler */
void atom::
zero_force ()
{
  f[0] = 0.0;
  f[1] = 0.0;
  f[2] = 0.0;
}

double atom::
kinetic_energy ()
{
  double vsq;
  
  vsq = v[0] * v[0] + v[1] * v[1] +
    v[2] * v[2];
  return 0.5 * mass * vsq;
}

static struct SPECIES_INFO
{
  enum SPECIES species;
  double mass, evdw, rvdw;	/* 10^-27 kg, maJ, 10^-10 m */
}
species_info[] =
{
  /*     10^-27 kg  10^-21 J  nm */
  C_SP3, 19.925, 0.357, 0.190,
  C_SP2, 19.925, 0.357, 0.194,
  C_SP, 19.925, 0.357, 0.194,
  H, 1.674, 0.382, 0.150,
  N_SP3, 23.251, 0.447, 0.182,
  O_SP3, 26.565, 0.406, 0.174,
  O_SP2, 26.565, 0.536, 0.174
};

static int num_species = sizeof (species_info) / sizeof (species_info[0]);

atom::atom (void)
{
  species = species_info[0].species;
  mass = species_info[0].mass;
  evdw = species_info[0].evdw;
  rvdw = species_info[0].rvdw;
  charge = 0.0;
  x[0] = x[1] = x[2] = 0.0;
  v[0] = v[1] = v[2] = 0.0;
  f[0] = f[1] = f[2] = 0.0;
}

atom::atom (enum SPECIES species, double charge,
	    double x0, double x1, double x2)
{
  int i, found;
  
  this->species = species;
  this->charge = charge;
  for (i = found = 0; i < num_species && !found; i++)
    if (species_info[i].species == species)
      {
	mass = species_info[i].mass;
	evdw = species_info[i].evdw;
	rvdw = species_info[i].rvdw;
	found = 1;
      }
  
  x[0] = x0;
  x[1] = x1;
  x[2] = x2;
  v[0] = v[1] = v[2] = 0.0;
  f[0] = f[1] = f[2] = 0.0;
}