mytypes.h

/*----------------------------------------------------------------------
  SerialReax - Reax Force Field Simulator

  Copyright (2010) Purdue University
  Hasan Metin Aktulga, haktulga@cs.purdue.edu
  Joseph Fogarty, jcfogart@mail.usf.edu
  Sagar Pandit, pandit@usf.edu
  Ananth Y Grama, ayg@cs.purdue.edu

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of
  the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  See the GNU General Public License for more details:
  <http://www.gnu.org/licenses/>.
  ----------------------------------------------------------------------*/

#ifndef __MYTYPES_H_
#define __MYTYPES_H_

#if (defined(HAVE_CONFIG_H) && !defined(__CONFIG_H_))
  #define __CONFIG_H_
  #include "config.h"
#endif

#include "math.h"
#include "random.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "sys/time.h"
#include "time.h"
#include "zlib.h"

#ifdef _OPENMP
  #include <omp.h>
#endif

//#define DEBUG_FOCUS
//#define TEST_FORCES
//#define TEST_ENERGY
#define REORDER_ATOMS  // turns on nbrgen opt by re-ordering atoms
//#define LGJ

#define SUCCESS  1
#define FAILURE  0
#define TRUE  1
#define FALSE 0

#define LOG    log
#define EXP    exp
#define SQRT   sqrt
#define POW    pow
#define ACOS   acos
#define COS    cos
#define SIN    sin
#define TAN    tan
#define CEIL   ceil
#define FLOOR  floor
#define FABS   fabs
#define FMOD   fmod

#define SQR(x)        ((x)*(x))
#define CUBE(x)       ((x)*(x)*(x))
#define DEG2RAD(a)    ((a)*PI/180.0)
#define RAD2DEG(a)    ((a)*180.0/PI)
#define MAX( x, y )   (((x) > (y)) ? (x) : (y))
#define MIN( x, y )   (((x) < (y)) ? (x) : (y))

/* NaN IEEE 754 representation for C99 in math.h
 * Note: function choice must match REAL typedef below */
#ifdef NAN
  #define IS_NAN_REAL(a) (isnan(a))
#else
  #warn "No support for NaN"
  #define NAN_REAL(a) (0)
#endif

#if !defined(PI)
  #define PI            3.14159265
#endif
#define C_ele          332.06371
//#define K_B         503.398008   // kcal/mol/K
#define K_B             0.831687   // amu A^2 / ps^2 / K
#define F_CONV          (1e6 / 48.88821291 / 48.88821291)   // --> amu A / ps^2
#define E_CONV          0.002391   // amu A^2 / ps^2 --> kcal/mol
#define EV_to_KCALpMOL 14.400000   // ElectronVolt --> KCAL per MOLe
#define KCALpMOL_to_EV 23.060549   // 23.020000//KCAL per MOLe --> ElectronVolt
#define ECxA_to_DEBYE   4.803204      // elem. charge * angstrom -> debye conv
#define CAL_to_JOULES   4.184000      // CALories --> JOULES
#define JOULES_to_CAL   (1/4.184000)    // JOULES --> CALories
#define AMU_to_GRAM     1.6605e-24
#define ANG_to_CM       1.0e-8
#define AVOGNR          6.0221367e23
#define P_CONV          (1.0e-24 * AVOGNR * JOULES_to_CAL)

#define MAX_STR             1024
#define MAX_LINE            1024
#define MAX_TOKENS          1024
#define MAX_TOKEN_LEN       1024

#define MAX_ATOM_ID         100000
#define MAX_RESTRICT        15
#define MAX_MOLECULE_SIZE   20
#define MAX_ATOM_TYPES      25

#define MAX_GRID            50
#define MAX_3BODY_PARAM     5
#define MAX_4BODY_PARAM     5
#define NO_OF_INTERACTIONS  10

#define MAX_dV              1.01
#define MIN_dV              0.99
#define MAX_dT              4.00
#define MIN_dT              0.00

#define ZERO           0.000000000000000e+00
#define ALMOST_ZERO    1e-10
#define NEG_INF       -1e10
#define NO_BOND        1e-3
#define HB_THRESHOLD   1e-2
#define MAX_BONDS      40
#define MIN_BONDS      15
#define MIN_HBONDS     50
#define SAFE_HBONDS    1.4
#define MIN_GCELL_POPL 50
#define SAFE_ZONE   1.2
#define DANGER_ZONE 0.95
#define LOOSE_ZONE  0.75


/* config params */
enum ensemble
{
    NVE = 0,
    NVT = 1,
    NPT = 2,
    sNPT = 3,
    iNPT = 4,
    ensNR = 5,
    bNVT = 6,
};

enum interaction_list_offets
{
    FAR_NBRS = 0,
    NEAR_NBRS = 1,
    THREE_BODIES = 2,
    BONDS = 3,
    OLD_BONDS = 4,
    HBONDS = 5,
    DBO = 6,
    DDELTA = 7,
    LIST_N = 8,
};

enum interaction_type
{
    TYP_VOID = 0,
    TYP_THREE_BODY = 1,
    TYP_BOND = 2,
    TYP_HBOND = 3,
    TYP_DBO = 4,
    TYP_DDELTA = 5,
    TYP_FAR_NEIGHBOR = 6,
    TYP_NEAR_NEIGHBOR = 7,
    TYP_N = 8,
};

enum errors
{
    FILE_NOT_FOUND = -10,
    UNKNOWN_ATOM_TYPE = -11,
    CANNOT_OPEN_FILE = -12,
    CANNOT_INITIALIZE = -13,
    INSUFFICIENT_MEMORY = -14,
    UNKNOWN_OPTION = -15,
    INVALID_INPUT = -16,
    INVALID_GEO = -17,
    NUMERIC_BREAKDOWN = -18,
    RUNTIME_ERROR = -19,
};

enum molecular_analysis_type
{
    NO_ANALYSIS = 0,
    FRAGMENTS = 1,
    REACTIONS = 2,
    NUM_ANALYSIS = 3,
};

enum restart_format
{
    WRITE_ASCII = 0,
    WRITE_BINARY = 1,
    RF_N = 2,
};

enum geo_formats
{
    CUSTOM = 0,
    PDB = 1,
    BGF = 2,
    ASCII_RESTART = 3,
    BINARY_RESTART = 4,
    GF_N = 5,
};

enum charge_method
{
    QEQ_CM = 0,
    EE_CM = 1,
    ACKS2_CM = 2,
};

enum solver
{
    GMRES_S = 0,
    GMRES_H_S = 1,
    CG_S = 2,
    SDM_S = 3,
};

enum pre_comp
{
    NONE_PC = 0,
    DIAG_PC = 1,
    ICHOLT_PC = 2,
    ILU_PAR_PC = 3,
    ILUT_PAR_PC = 4,
    ILU_SUPERLU_MT_PC = 5,
    SAI_PC = 6,
};

enum pre_app
{
    TRI_SOLVE_PA = 0,
    TRI_SOLVE_LEVEL_SCHED_PA = 1,
    TRI_SOLVE_GC_PA = 2,
    JACOBI_ITER_PA = 3,
};


typedef double real;
typedef real rvec[3];
typedef int ivec[3];
typedef real rtensor[3][3];


/* Force field global params mapping:
 *
 * l[0]  = p_boc1
 * l[1]  = p_boc2
 * l[2]  = p_coa2
 * l[3]  = N/A
 * l[4]  = N/A
 * l[5]  = N/A
 * l[6]  = p_ovun6
 * l[7]  = N/A
 * l[8]  = p_ovun7
 * l[9]  = p_ovun8
 * l[10] = N/A
 * l[11] = N/A
 * l[12] = N/A
 * l[13] = N/A
 * l[14] = p_val6
 * l[15] = p_lp1
 * l[16] = p_val9
 * l[17] = p_val10
 * l[18] = N/A
 * l[19] = p_pen2
 * l[20] = p_pen3
 * l[21] = p_pen4
 * l[22] = N/A
 * l[23] = p_tor2
 * l[24] = p_tor3
 * l[25] = p_tor4
 * l[26] = N/A
 * l[27] = p_cot2
 * l[28] = p_vdW1
 * l[29] = v_par30
 * l[30] = p_coa4
 * l[31] = p_ovun4
 * l[32] = p_ovun3
 * l[33] = p_val8
 * l[34] = ACKS2 bond softness
 * l[35] = N/A
 * l[36] = N/A
 * l[37] = version number
 * l[38] = p_coa3 */
typedef struct
{
    int n_global;
    real* l;
    int vdw_type;
} global_parameters;


typedef struct
{
    /* Line one in field file */
    char name[15];                     /* Two character atom name */

    real r_s;
    real valency;                     /* Valency of the atom */
    real mass;                        /* Mass of atom */
    real r_vdw;
    real epsilon;
    real gamma;
    real r_pi;
    real valency_e;
    real nlp_opt;

    /* Line two in field file */
    real alpha;
    real gamma_w;
    real valency_boc;
    real p_ovun5;
    real chi;
    real eta;
    /* Determines whether this type of atom participates in H_bonds.
     * It is 1 for donor H, 2 for acceptors (O,S,N), 0 for others*/
    int p_hbond;

    /* Line three in field file */
    real r_pi_pi;
    real p_lp2;
    real b_o_131;
    real b_o_132;
    real b_o_133;
    /* bond softness for ACKS2 */
    real b_s_acks2;

    /* Line four in the field file */
    real p_ovun2;
    real p_val3;
    real valency_val;
    real p_val5;
    real rcore2;
    real ecore2;
    real acore2;
} single_body_parameters;


/* Two Body Parameters */
typedef struct
{
    /* Bond Order parameters */
    real p_bo1;
    real p_bo2;
    real p_bo3;
    real p_bo4;
    real p_bo5;
    real p_bo6;
    real r_s;
    real r_p;
    real r_pp;  /* r_o distances in BO formula */
    real p_boc3;
    real p_boc4;
    real p_boc5;

    /* Bond Energy parameters */
    real p_be1;
    real p_be2;
    real De_s;
    real De_p;
    real De_pp;

    /* Over/Under coordination parameters */
    real p_ovun1;

    /* Van der Waal interaction parameters */
    real D;
    real alpha;
    real r_vdW;
    real gamma_w;
    real rcore;
    real ecore;
    real acore;

    /* electrostatic parameters */
    real gamma; // note: this parameter is gamma^-3 and not gamma.

    real v13cor, ovc;
} two_body_parameters;


/* 3-body parameters */
typedef struct
{
    /* valence angle */
    real theta_00;
    real p_val1;
    real p_val2;
    real p_val4;
    real p_val7;

    /* penalty */
    real p_pen1;

    /* 3-body conjugation */
    real p_coa1;
} three_body_parameters;


typedef struct
{
    int cnt;
    three_body_parameters prm[MAX_3BODY_PARAM];
} three_body_header;


/* hydrogen-bond parameters */
typedef struct
{
    real r0_hb;
    real p_hb1;
    real p_hb2;
    real p_hb3;
} hbond_parameters;


/* 4-body parameters */
typedef struct
{
    real V1;
    real V2;
    real V3;

    /* torsion angle */
    real p_tor1;

    /* 4-body conjugation */
    real p_cot1;
} four_body_parameters;


typedef struct
{
    int cnt;
    four_body_parameters prm[MAX_4BODY_PARAM];
} four_body_header;


typedef struct
{
    int num_atom_types;
    global_parameters gp;
    single_body_parameters *sbp;
    two_body_parameters **tbp;
    three_body_header ***thbp;
    hbond_parameters ***hbp;
    four_body_header ****fbp;
} reax_interaction;


typedef struct
{
    /* Type of this atom */
    int type;
    /**/
    char name[8];
    /* position */
    rvec x;
    /* velocity */
    rvec v;
    /* force */
    rvec f;
    /* Charge on the atom */
    real q;
} reax_atom;


typedef struct
{
    real volume;

    rvec box_norms;
    rvec side_prop;
    rvec nbr_box_press[27];

    rtensor box, box_inv, old_box;
    rtensor trans, trans_inv;
    rtensor g;
} simulation_box;


typedef struct
{
    int max_atoms;
    int max_nbrs;
    int total;
    real cell_size;
    ivec spread;

    ivec ncell;
    rvec len;
    rvec inv_len;

    int**** atoms;
    int*** top;
    int*** mark;
    int*** start;
    int*** end;
    ivec**** nbrs;
    rvec**** nbrs_cp;
} grid;


typedef struct
{
    /* number of atoms */
    int N;
    /* dimension of the N x N sparse charge method matrix H */
    int N_cm;
    /* atom info */
    reax_atom *atoms;
    /* atomic interaction parameters */
    reax_interaction reaxprm;
    /* simulation space (a.k.a. box) parameters */
    simulation_box box;
    /* grid structure used for binning atoms and tracking neighboring bins */
    grid g;
} reax_system;


/* Simulation control parameters not related to the system */
typedef struct
{
    char sim_name[MAX_STR];
    char restart_from[MAX_STR];
    int restart;
    int random_vel;

    int reposition_atoms;

    /* ensemble values:
       0 : NVE
       1 : NVT  (Nose-Hoover)
       2 : NPT  (Parrinello-Rehman-Nose-Hoover) Anisotropic
       3 : sNPT (Parrinello-Rehman-Nose-Hoover) semiisotropic
       4 : iNPT (Parrinello-Rehman-Nose-Hoover) isotropic */
    int ensemble;
    int nsteps;
    int periodic_boundaries;
    int restrict_bonds;
    int tabulate;
    ivec periodic_images;
    real dt;

    int reneighbor;
    real vlist_cut;
    real nbr_cut;
    real r_cut;
    real r_sp_cut;
    real r_low; // upper and lower taper
    real bo_cut;
    real thb_cut;
    real hb_cut;
    real Tap7;
    real Tap6;
    real Tap5;
    real Tap4;
    real Tap3;
    real Tap2;
    real Tap1;
    real Tap0;

    real T_init;
    real T_final;
    real T;
    real Tau_T;
    int T_mode;
    real T_rate;
    real T_freq;

    real Tau_PT;
    rvec P;
    rvec Tau_P;
    int press_mode;
    real compressibility;

    int remove_CoM_vel;

    int geo_format;

    int dipole_anal;
    int freq_dipole_anal;

    int diffusion_coef;
    int freq_diffusion_coef;
    int restrict_type;

    unsigned int charge_method;
    unsigned int cm_solver_type;
    real cm_q_net;
    unsigned int cm_solver_max_iters;
    unsigned int cm_solver_restart;
    real cm_solver_q_err;
    real cm_domain_sparsity;
    unsigned int cm_domain_sparsify_enabled;
    unsigned int cm_solver_pre_comp_type;
    unsigned int cm_solver_pre_comp_refactor;
    real cm_solver_pre_comp_droptol;
    unsigned int cm_solver_pre_comp_sweeps;
    real cm_solver_pre_comp_sai_thres;
    unsigned int cm_solver_pre_app_type;
    unsigned int cm_solver_pre_app_jacobi_iters;

    int molec_anal;
    int freq_molec_anal;
    real bg_cut;
    int num_ignored;
    int ignore[MAX_ATOM_TYPES];

    int num_threads;
} control_params;


typedef struct
{
    real T;
    real xi;
    real v_xi;
    real v_xi_old;
    real G_xi;
} thermostat;


typedef struct
{
    real P;
    real eps;
    real v_eps;
    real v_eps_old;
    real a_eps;

} isotropic_barostat;


typedef struct
{
    rtensor P;
    real P_scalar;

    real eps;
    real v_eps;
    real v_eps_old;
    real a_eps;

    rtensor h0;
    rtensor v_g0;
    rtensor v_g0_old;
    rtensor a_g0;

} flexible_barostat;


typedef struct
{
    real start;
    real end;
    real elapsed;

    real total;
    real nbrs;
    real init_forces;
    real bonded;
    real nonb;
    real cm;
    real cm_sort_mat_rows;
    real cm_solver_pre_comp;
    real cm_solver_pre_app;
    int cm_solver_iters;
    real cm_solver_spmv;
    real cm_solver_vector_ops;
    real cm_solver_orthog;
    real cm_solver_tri_solve;
} reax_timing;


typedef struct
{
    int step;
    int prev_steps;
    real time;

    real M;              /* Total Mass */
    real inv_M;                      /* 1 / Total Mass */

    rvec xcm;                        /* Center of mass */
    rvec vcm;                        /* Center of mass velocity */
    rvec fcm;                        /* Center of mass force */
    rvec amcm;                       /* Angular momentum of CoM */
    rvec avcm;                       /* Angular velocity of CoM */
    real etran_cm;                  /* Translational kinetic energy of CoM */
    real erot_cm;                    /* Rotational kinetic energy of CoM */

    rtensor kinetic;                 /* Kinetic energy tensor */
    rtensor virial;                  /* Hydrodynamic virial */

    real E_Tot;
    real E_Kin;                      /* Total kinetic energy */
    real E_Pot;

    real E_BE;                       /* Total bond energy */
    real E_Ov;                       /* Total over coordination */
    real E_Un;                       /* Total under coordination energy */
    real E_Lp;                       /* Total under coordination energy */
    real E_Ang;                      /* Total valance angle energy */
    real E_Pen;                      /* Total penalty energy */
    real E_Coa;                      /* Total three body conjgation energy */
    real E_HB;                       /* Total Hydrogen bond energy */
    real E_Tor;                      /* Total torsional energy */
    real E_Con;                      /* Total four body conjugation energy */
    real E_vdW;                      /* Total van der Waals energy */
    real E_Ele;                      /* Total electrostatics energy */
    real E_Pol;                      /* Polarization energy */

    real N_f;                        /*Number of degrees of freedom */
    rvec t_scale;
    rtensor p_scale;
    thermostat therm;                /* Used in Nose_Hoover method */
    isotropic_barostat iso_bar;
    flexible_barostat flex_bar;
    real inv_W;

    rvec int_press;
    rvec ext_press;
    real kin_press;
    rvec tot_press;

    reax_timing timing;
} simulation_data;


typedef struct
{
    int thb;
    /* pointer to the third body on the central atom's nbrlist */
    int pthb;
    real theta;
    real cos_theta;
    rvec dcos_di;
    rvec dcos_dj;
    rvec dcos_dk;
} three_body_interaction_data;


typedef struct
{
    int nbr;
    ivec rel_box;
//    rvec ext_factor;
    real d;
    rvec dvec;
} near_neighbor_data;


typedef struct
{
    int nbr;
    ivec rel_box;
//    rvec ext_factor;
    real d;
    rvec dvec;
} far_neighbor_data;


typedef struct
{
    int nbr;
    int scl;
    far_neighbor_data *ptr;
} hbond_data;


typedef struct
{
    int wrt;
    rvec dVal;
} dDelta_data;


typedef struct
{
    int wrt;
    rvec dBO;
    rvec dBOpi;
    rvec dBOpi2;
} dbond_data;


typedef struct
{
    real BO;
    real BO_s;
    real BO_pi;
    real BO_pi2;
    real Cdbo;
    real Cdbopi;
    real Cdbopi2;
    real C1dbo;
    real C2dbo;
    real C3dbo;
    real C1dbopi;
    real C2dbopi;
    real C3dbopi;
    real C4dbopi;
    real C1dbopi2;
    real C2dbopi2;
    real C3dbopi2;
    real C4dbopi2;
    rvec dBOp;
    rvec dln_BOp_s;
    rvec dln_BOp_pi;
    rvec dln_BOp_pi2;
} bond_order_data;


typedef struct
{
    int nbr;
    int sym_index;
    int dbond_index;
    ivec rel_box;
//    rvec ext_factor;
    real d;
    rvec dvec;
    bond_order_data bo_data;
} bond_data;


/* compressed row storage (crs) format
 * See, e.g.,
 *   http://netlib.org/linalg/html_templates/node91.html#SECTION00931100000000000000
 *
 *   m: number of nonzeros (NNZ) ALLOCATED
 *   n: number of rows
 *   start: row pointer (last element contains ACTUAL NNZ)
 *   j: column index for corresponding matrix entry
 *   val: matrix entry */
typedef struct
{
    unsigned int n;
    unsigned int m;
    unsigned int *start;
    unsigned int *j;
    real *val;
} sparse_matrix;


typedef struct
{
    int num_far;
    int Htop;
    int hbonds;
    int num_hbonds;
    int bonds;
    int num_bonds;
    int num_3body;
    int gcell_atoms;
} reallocate_data;


typedef struct
{
    real H;
    real e_vdW;
    real CEvd;
    real e_ele;
    real CEclmb;
} LR_data;


typedef struct
{
    real a;
    real b;
    real c;
    real d;
} cubic_spline_coef;


typedef struct
{
    real xmin;
    real xmax;
    int n;
    real dx;
    real inv_dx;
    real a;

    real m;
    real c;

    real *y;
} lookup_table;


typedef struct
{
    real xmin;
    real xmax;
    int n;
    real dx;
    real inv_dx;
    real a;
    real m;
    real c;

    LR_data *y;
    cubic_spline_coef *H;
    cubic_spline_coef *vdW;
    cubic_spline_coef *CEvd;
    cubic_spline_coef *ele;
    cubic_spline_coef *CEclmb;
} LR_lookup_table;


typedef struct
{
    /* bond order related storage */
    real *total_bond_order;
    real *Deltap;
    real *Deltap_boc;
    real *Delta;
    real *Delta_lp;
    real *Delta_lp_temp;
    real *Delta_e;
    real *Delta_boc;
    real *dDelta_lp;
    real *dDelta_lp_temp;
    real *nlp;
    real *nlp_temp;
    real *Clp;
    real *vlpex;
    rvec *dDeltap_self;

    /* charge method storage */
    sparse_matrix *H;
    sparse_matrix *H_full;
    sparse_matrix *H_sp;
    sparse_matrix *H_p;
    sparse_matrix *H_spar_patt;
    sparse_matrix *H_spar_patt_full;
    sparse_matrix *H_app_inv;
    sparse_matrix *L;
    sparse_matrix *U;
    real *droptol;
    real *Hdia_inv;
    real *b;
    real *b_s;
    real *b_t;
    real *b_prc;
    real *b_prm;
    real **s;
    real **t;

    /* GMRES related storage */
    real *y;
    real *z;
    real *g;
    real *hc;
    real *hs;
    real **h;
    real **rn;
    real **v;

    /* CG related storage */
    real *r;
    real *d;
    real *q;
    real *p;

    /* SpMV related storage */
#ifdef _OPENMP
    real *b_local;
#endif

    /* Level scheduling related storage for applying, e.g. ICHOLT and ILU(T),
     * preconditioners */
    int levels_L;
    int levels_U;
    unsigned int *row_levels_L;
    unsigned int *level_rows_L;
    unsigned int *level_rows_cnt_L;
    unsigned int *row_levels_U;
    unsigned int *level_rows_U;
    unsigned int *level_rows_cnt_U;
    unsigned int *top;

    /* Graph coloring related storage for applying, e.g. ICHOLT and ILU(T),
     * preconditioners */
    unsigned int *color;
    unsigned int *to_color;