#ifndef _MD_H
#define _MD_H
#include <NIST.h>
#include <Min.h>
#include <Atoms.h>
#include <Timer.h>
#include <Phonon.h>
#include <SSpectral.h>

#define Zr_ATOM        0
#define C_ATOM         1
#define MAX_TYPE       2
/* #define RCUT_ZrZr      7.0 */
#define RCUT_ZrZr      6.1
#define K_ZrZr         0.1
#define RCUT_EDGE     (K_ZrZr / 100)
#define RCUT_ZrZr_0   (RCUT_ZrZr-RCUT_EDGE)
#define A_ZrZr         2.9296875
#define B_ZrZr         2.58787395638939
#define C_ZrZr         1.32308467741935
#define D_ZrZr         4.3672464262563
#define RCUT_ZrC       3.5
#define RCUT_ZrC_0    (RCUT_ZrC-RCUT_EDGE)
#define K_ZrC          K_ZrZr
#define A_ZrC          3.24589393669854
#define B_ZrC          2.05679804919117
#define C_ZrC          0.82303818052368
#define D_ZrC          4.15482225815134
#define LAMBDA         0
#define ALPHA_ZrC      1.80853303846249
#define BETA_ZrC       14.59345494373451
#define ALPHA_CZr      ALPHA_ZrC
#define BETA_CZr       BETA_ZrC
#define RCUT_CC        0
#define MAX_FIRST      16
#define MAX_FIRST4     64

#define DEF_NEIGHBORMAX             128         /* before pairwise-saving */
#define DEF_MIN_SHRINKAGE           0.99        /* compared to session H0 */
#define DEF_MAX_ATOM_DISPLACEMENT   0.35        /* Lindemann's law */
#define DEF_TMP_CONFIG_FREQ         20          /* x report frequency */
#define DEF_COUPLING_TIME_CONSTANT  300         /* x g */
#define FRESHSTART_OVERDRIVE        2.0
#define JUMPSTART_T_IN_K            0.1         /* initialize a spark */
#define NTS_MAX_WALL_KINE_IN_KT     15.         /* prevent explosion */
#define NTS_CALM_AND_START_AVERAGE  0.5
#define NTS_CALM_AND_STOP_AVERAGE   0.95

typedef struct
{
    double real;
    double md;
    double dTmd__dTreal;
} Temperature;

typedef struct
{
    double counter;
    M3 H;
} Geometry;

typedef struct
{
    int NTS_calm_in_g;
    int NTS_calm_and_start_average_in_g;
    int NTS_calm_and_stop_average_in_g;
    bool do_preliminary;
    bool preliminary_done;
    Geometry geometry;
} Preliminary;

typedef struct
{
    bool snapshot_taken;
    int t_in_g;
    double MSD;
} Snapshot;

typedef struct
{
    double counter;
    M3 H;
    double kine;
    double pote;
    M3 stress;
} Property;

#define drawtime (t_in_g * g * UTIME_IN_S)
#define draw(name,value) \
  fprintf( fp_draw, name" %e %g\n", drawtime, (double)(value) )
#define DRAW(NAME,value) draw( STR(NAME), value )
#define M3DRAW(NAME,A) (       draw(STR(NAME)"11",A[0][0]), \
  draw(STR(NAME)"12",A[0][1]), draw(STR(NAME)"13",A[0][2]), \
  draw(STR(NAME)"21",A[1][0]), draw(STR(NAME)"22",A[1][1]), \
  draw(STR(NAME)"23",A[1][2]), draw(STR(NAME)"31",A[2][0]), \
  draw(STR(NAME)"32",A[2][1]), draw(STR(NAME)"33",A[2][2]) )
#define M3DRAWMUL(NAME,A,factor) (      draw(STR(NAME)"11",A[0][0]*(factor)),\
  draw(STR(NAME)"12",A[0][1]*(factor)), draw(STR(NAME)"13",A[0][2]*(factor)),\
  draw(STR(NAME)"21",A[1][0]*(factor)), draw(STR(NAME)"22",A[1][1]*(factor)),\
  draw(STR(NAME)"23",A[1][2]*(factor)), draw(STR(NAME)"31",A[2][0]*(factor)),\
  draw(STR(NAME)"32",A[2][1]*(factor)), draw(STR(NAME)"33",A[2][2]*(factor)) )
#define SYMMAT_DRAW(NAME,A) ( \
  draw(STR(NAME)"11",A[0][0]), draw(STR(NAME)"12",A[0][1]), \
  draw(STR(NAME)"13",A[0][2]), draw(STR(NAME)"22",A[1][1]), \
  draw(STR(NAME)"23",A[1][2]), draw(STR(NAME)"33",A[2][2]) )
#define SYMMAT_DRAWMUL(NAME,A,factor) ( \
  draw(STR(NAME)"11",A[0][0]*(factor)), draw(STR(NAME)"12",A[0][1]*(factor)),\
  draw(STR(NAME)"13",A[0][2]*(factor)), draw(STR(NAME)"22",A[1][1]*(factor)),\
  draw(STR(NAME)"23",A[1][2]*(factor)), draw(STR(NAME)"33",A[2][2]*(factor)) )

#endif  /* _MD_H */