hawaii.h

Go to the documentation of this file.

#ifndef HAWAII_H_
#define HAWAII_H_
 
#include <assert.h>
#include <errno.h>
#include <stdarg.h>
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <regex.h>
#include <time.h>
#include <unistd.h>
 
#ifdef USE_MPI
#include <mpi.h>
#endif // USE_MPI
 
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_spline.h>
 
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>
 
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
 
#define OMPI_SKIP_MPICXX
 
#ifndef __cplusplus
#define _GNU_SOURCE
void sincos(double, double*, double*);
int syncfs(int);
#endif
#include <math.h>
#include <complex.h>
 
#ifndef __cplusplus
#define MT_GENERATE_CODE_IN_HEADER 0
#endif 
#include "mtwist.h"
/*
 * We are using the following functions from mtwist:
 * double mt_drand(void)
 *   Return a pseudorandom double in [0,1) with 32 bits of randomness
 *
 * uint32_t mt_lrand(void);
 *   Generate 32-bit random value 
 */
 
 
#include <cvode/cvode.h>
#include <nvector/nvector_serial.h>
 
#include "array.h"
#include "constants.h"
#include "arena.h"
 
#define IPHI    0 
#define IPPHI   1 
#define ITHETA  2 
#define IPTHETA 3 
#define IR      4 
#define IPR     5 
// used only in Monomer.qp
#define IPSI    4 
#define IPPSI   5 
 
#define UNUSED(x) (void)(x)
#define TODO(message) do { fprintf(stderr, "%s:%d: TODO: %s\n", __FILE__, __LINE__, message); abort(); } while(0)
#define UNREACHABLE(message) do { fprintf(stderr, "%s:%d: UNREACHABLE: %s\n", __FILE__, __LINE__, message); abort(); } while(0)
 
#define return_defer(value) do { result = (value); goto defer; } while(0)
 
#define DA_INIT_CAP 256
#define da_append(da, item)                                                          \
    do {                                                                             \
        if ((da)->count >= (da)->capacity) {                                         \
            (da)->capacity = (da)->capacity == 0 ? DA_INIT_CAP : (da)->capacity*2;   \
            (da)->items = realloc((da)->items, (da)->capacity*sizeof(*(da)->items)); \
            assert((da)->items != NULL && "ASSERT: not enough memory\n");            \
        }                                                                            \
                                                                                     \
        (da)->items[(da)->count++] = (item);                                         \
    } while (0)
 
#define da_insert(da, i, item) \
    do {                                 \
        if ((i < 0) || ((i) > (da)->count)) { \
            assert(0 && "ASSERT: index out of bounds\n"); \
        } \
        if ((da)->count >= (da)->capacity) { \
            (da)->capacity = (da)->capacity == 0 ? DA_INIT_CAP : (da)->capacity*2; \
            (da)->items = realloc((da)->items, (da)->capacity*sizeof(*(da)->items)); \
            assert((da)->items != NULL && "ASSERT: not enough memory\n"); \
        } \
        memmove((da)->items + (i) + 1, (da)->items + (i), ((da)->count - (i))*sizeof(*(da)->items)); \
        (da)->items[(i)] = (item); \
        (da)->count++; \
    } while(0)
 
#define da_last(da) (da)->items[(da)->count - 1]
 
extern int  _wrank;
extern int  _wsize;
extern bool _print0_suppress_info;
extern int  _print0_margin; 
 
#ifdef USE_MPI
#define INIT_WRANK                          \
    MPI_Comm_rank(MPI_COMM_WORLD, &_wrank); \
 
#define INIT_WSIZE                          \
    MPI_Comm_size(MPI_COMM_WORLD, &_wsize); \
 
#define INFO(...)                                                   \
    if ((_wrank == 0) && !_print0_suppress_info) {                  \
        if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
        printf("INFO: "); printf(__VA_ARGS__);                      \
    }
 
#define WARNING(...)                                                \
    if (_wrank == 0) {                                              \
        if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
        printf("WARNING: "); printf(__VA_ARGS__);                   \
    }
 
#define ERROR(...) \
    if (_wrank == 0) {                                            \
      if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
      printf("ERROR: "); printf(__VA_ARGS__);                     \
    } 
 
#define PRINT0(...)                                               \
    if (_wrank == 0) {                                            \
      if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
      printf(__VA_ARGS__);                                        \
    }
 
#else
#define INIT_WRANK
#define INIT_WSIZE
 
#define INFO(...)                                                   \
    if (!_print0_suppress_info) {                                   \
        if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
        printf(__VA_ARGS__);                                        \
    }
 
#define WARNING(...)                                            \
    if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
    printf("WARNING: "); printf(__VA_ARGS__);                   \
 
#define ERROR(...)                                              \
    if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
    printf("ERROR: "); printf(__VA_ARGS__);                     \
 
#define PRINT0(...)                                             \
    if (_print0_margin > 0) printf("%*s", _print0_margin, " "); \
    printf(__VA_ARGS__);
 
 
#endif
 
#ifdef __cplusplus
extern "C" {
#endif
 
#define MONOMER_COUNT 6 
#define MODULO_BASE 100
typedef enum {
    ATOM                                 = 0,                 
    LINEAR_MOLECULE                      = 4,                 
    LINEAR_MOLECULE_REQ_INTEGER          = MODULO_BASE + 4,   
    LINEAR_MOLECULE_REQ_HALFINTEGER      = 2*MODULO_BASE + 4, 
    ROTOR                                = 6,                 
    ROTOR_REQUANTIZED_ROTATION           = MODULO_BASE + 6,   
} MonomerType;
 
extern MonomerType MONOMER_TYPES[MONOMER_COUNT];
 
// note: gsl_histogram does have fields for nbins and max but I don't want
// to access them when deciding on what the values should be, so we will have "duplicates"
// The values of 'nbins' and 'max' are the initial values, and can be updated when 
// the histogram is extended. 
typedef struct {
    gsl_histogram *h;
    size_t nbins;  
    double max;
    char *filename;
    FILE *fp;
    bool is_allocated;
} StoredHistogram;

typedef struct { 
    int index;                 
    MonomerType t;             
    double II[3];              
    double DJ;                 
    
    double *qp;                
    double *dVdq;              
   
    bool apply_requantization; 
   
    size_t req_switch_counter; 
    size_t torque_cache_len;   
    double torque_limit;       
    double *torque_cache;      
 
    double initial_j;
    
    /* histograms */
    StoredHistogram nswitch_histogram;
 
    size_t jini_histogram_bins;
    double jini_histogram_max;
    char *jini_histogram_filename; 
    gsl_histogram *jini_histogram;
    FILE *fp_jini_histogram; 
   
    size_t jfin_histogram_bins;
    double jfin_histogram_max;
    char *jfin_histogram_filename;
    gsl_histogram *jfin_histogram;
    FILE *fp_jfin_histogram;
} Monomer; 

typedef struct {
    double intermolecular_qp[6]; 
    Monomer m1;                  
    Monomer m2;                  
    double mu;                   
 
    size_t Q_SIZE;               
    size_t QP_SIZE;              
 
    double *intermediate_q;      
    double *dVdq;                
 
    size_t seed;
 
    // time_t: time as the number of seconds since the Epoch (1970-01-01)  
    time_t init_rawtime;         
    time_t temp_rawtime;         
} MoleculeSystem;

typedef enum {
    PAIR_STATE_NONE,                
    PAIR_STATE_FREE_AND_METASTABLE, 
    PAIR_STATE_BOUND,               
    PAIR_STATE_ALL,                 
    PAIR_STATE_COUNT,               
} PairState;
 
extern const char *PAIR_STATES[PAIR_STATE_COUNT]; 
 
typedef enum {
    CALCULATION_NONE,
    CALCULATION_PR_MU,
    CALCULATION_CORRELATION_SINGLE,
    CALCULATION_CORRELATION_ARRAY,
    CALCULATION_PROCESSING,
    CALCULATION_PHASE_SPACE_M0,
    CALCULATION_PHASE_SPACE_M2,
    CALCULATION_TYPES_COUNT,
} CalculationType;
 
extern const char* CALCULATION_TYPES[CALCULATION_TYPES_COUNT];
 
typedef struct {
    PairState ps;
    CalculationType calculation_type;
 
    /* sampling */
    double sampler_Rmin; // a.u.
    double sampler_Rmax; // a.u.
    double pesmin; // Hartree 
 
    /* initial spectral moments check */ 
    size_t initialM0_npoints;
    size_t initialM2_npoints;
    double partial_partition_function_ratio;
 
    size_t hep_m0_niterations;
    size_t hep_m0_npoints;
    size_t hep_m2_niterations;
    size_t hep_m2_npoints;
    size_t hep_ppf_niterations;
    size_t hep_ppf_npoints;
 
    /* weights to factor in spin statistics */
    double odd_j_spin_weight;
    double even_j_spin_weight;
 
    /* trajectory */
    double sampling_time; // a.t.u.
    size_t MaxTrajectoryLength;
    bool allow_truncating_trajectories_at_length_limit;
    double cvode_tolerance;
   
    /* iteration parameters :: applicable to both correlation function AND spectral function calculations */ 
    size_t niterations; // the accumulation of the total_trajectories is split into this number of iterations
    size_t total_trajectories; // the total number of accumulated trajectories (ALL iterations) 
 
    /* correlation function and correlation function array calculations ONLY */
    const char* cf_filename;
    double Rcut; // distance at which the trajectory is forcefully stopped, a.u. 
    bool use_zimmermann_trick;
 
    /* pr/mu spectral function calculation ONLY */   
    const char *sf_filename;
    double ApproximateFrequencyMax; // cm-1
    double R0; // initial distance, a.u. 
    double average_time_between_collisions; // a.t.u.
 
    /* correlation function array ONLY */
    double* partial_partition_function_ratios;
    double *satellite_temperatures;
    size_t num_satellite_temperatures;
    const char **cf_filenames;
} CalcParams;
 
 
 
typedef enum {
    REJECT,
    MCMC
} SamplingType;
 
 
typedef struct {
    double *t; 
    double *data;
    size_t len;      // # of samples in *t, *data
    size_t capacity; // capacity *t, *data
    double ntraj;    // # of trajectories used for averaging
    double Temperature;
    bool normalized; // flag that indicates whether the *data samples are normalized by # of trajectories 
} CFnc;
 
typedef struct {
    CFnc **items;
    size_t count;
    size_t capacity;
} CFncs;
 
 
/*
 * An array of correlation functions for a fixed 'base temperature'
 * 'base temperature'      -- a temperature for which sampling of initial conditions is performed
 * 'satellite temperature' -- a target temperature for which correlation function is re-weighted 
 */ 
typedef struct {
    double *t;
    double **data;
    size_t ntemp;
    size_t len;    // # of samples in *t and elements of *data 
    double *nstar; // array of # effective trajectories
    size_t ntraj; // # of calculated trajectories
} CFncArray;
 
// SFnc may originate from CFnc which itself was derived 
// via multi-temperature reweighting. As a result its 'ntraj' 
// turns out be floating-point value rather than an integer
typedef struct {
    double *nu;
    double *data;
    size_t len;      // # of samples in *nu, *data
    size_t capacity; // capacity of *nu, *data
    double ntraj;    // # of trajectories used for averaging
    double Temperature;
    bool normalized; 
} SFnc;
 
 
typedef struct {
    double *nu;
    double *data;
    size_t len;      // # of samples of *nu, *data
    size_t capacity; // capacity of *nu, *data
    double ntraj;    // # of trajectories used for averaging
    double Temperature;
    bool normalized;
} Spectrum;
 
typedef struct {
    Spectrum **items;
    size_t count;
    size_t capacity;
} Spectra;
 
/*
 * MODEL: Lorentzian function shifted upwards by constant: 
 *             y = C + A /(1 + B^2 x^2)
 */
typedef struct {
    double A;
    double B;
    double C;
} WingParams;
    
typedef struct {
    size_t n;
    double* t;
    double* y;
} WingData;
 
typedef struct {
    double before2;
    double before;
    double current;
 
    size_t turning_points;
    
    size_t called;
    bool ready;
} Tracker;
 
 
/* ---------------------------------------------------------------------------------------------------------------- */
/* --------------------------   User-Supplied Functions: loaded from dynamic libs --------------------------------- */
/* ---------------------------------------------------------------------------------------------------------------- */ 
typedef void (*dipolePtr)(double*, double[3]);
extern dipolePtr dipole_1;
extern dipolePtr dipole_2;
 
typedef void (*dipoleFree)(void);
extern dipoleFree free_dipole_1;
extern dipoleFree free_dipole_2;

typedef double (*pesPtr)(double*);
extern pesPtr pes;
 
typedef void (*dpesPtr)(double*, double*);
extern dpesPtr dpes;
/* ---------------------------------------------------------------------------------------------------------------- */
 
MoleculeSystem init_ms(double mu, MonomerType t1, MonomerType t2, double *I1, double *I2, size_t seed); 
MoleculeSystem init_ms_from_monomers(double mu, Monomer *m1, Monomer *m2, size_t seed);
void free_ms(MoleculeSystem *ms);
 
#ifdef USE_MPI
CFnc calculate_correlation_and_save(MoleculeSystem *ms, CalcParams *params, double Temperature);
SFnc calculate_spectral_function_using_prmu_representation_and_save(MoleculeSystem *ms, CalcParams *params, double Temperature);
CFncArray calculate_correlation_array_and_save(MoleculeSystem *ms, CalcParams *params, double base_temperature); 
#endif // USE_MPI
 
const char* display_monomer_type(MonomerType t);
 
void put_qp_into_ms(MoleculeSystem *ms, Array qp);
void get_qp_from_ms(MoleculeSystem *ms, Array *qp);
 
// this function takes the phase point from MoleculeSystem 
// and packs it into "intermediate_q" field. Then this pointer
// can be passed to potential energy / its derivatives / dipole function 
void extract_q_and_write_into_ms(MoleculeSystem *ms);
 
// this function takes the corresponding components of "ms.dVdq" and writes
// them into "m.dVdq" vectors for each monomer 
void extract_dVdq_and_write_into_monomers(MoleculeSystem *ms);
 
// 20.12.2024 NOTE: 
// the MoleculeSystem struct needs to be passed as void* into "rhs" function.
// Then, the first step is to write the phase point (N_Vector y) into the fields of MoleculeSystem
// using the method "put_qp_into_ms".
int rhs(realtype t, N_Vector y, N_Vector ydot, void *data);
 
Array compute_numerical_derivatives(double (*f)(double *q), double *q, size_t len, size_t order);
Array compute_numerical_rhs(MoleculeSystem *ms, size_t order); 
gsl_matrix* compute_numerical_jac(void (*transform_angles)(double *qlab, double *qmol), double *qlab, size_t ninput_coordinates, size_t noutput_coordinates, size_t order);
 
double kinetic_energy(MoleculeSystem *ms);
double Hamiltonian(MoleculeSystem *ms);
 
double generate_normal(double sigma); 
 
void q_generator(MoleculeSystem *ms, CalcParams *params);
void p_generator(MoleculeSystem *ms, double Temperature);
bool reject(MoleculeSystem *ms, double Temperature, double pesmin);
 
void try_applying_requantization_for_monomer(Monomer *m, size_t step_counter);
void j_monomer(Monomer *m, double j[3]);
double torque_monomer(Monomer *m); 
double find_closest_integer(double j);
double find_closest_half_integer(double j); 
 
void invert_momenta(MoleculeSystem *ms);
 
double analytic_full_partition_function_by_V(MoleculeSystem *ms, double Temperature);
 
// ----------------------------------------------------------
// Spectral moments calculation 
// ----------------------------------------------------------
double integrate_composite_simpson(double *x, double *y, size_t len); 
bool compute_Mn_from_cf_using_classical_detailed_balance(CFnc cf, size_t n, double *result);
double compute_Mn_from_sf_using_classical_detailed_balance(SFnc sf, size_t n);
double compute_Mn_from_sf_using_quantum_detailed_balance(SFnc sf, size_t n);
void calculate_M0(MoleculeSystem *ms, CalcParams *params, double Temperature, double *m, double *q);
void compute_dHdp(MoleculeSystem *ms, gsl_matrix* dHdp); 
void calculate_M2(MoleculeSystem *ms, CalcParams *params, double Temperature, double *m, double *q);
 
#ifdef USE_MPI
void mpi_calculate_M0(MoleculeSystem *ms, CalcParams *params, double Temperature, double *m, double *q);
void mpi_calculate_M2(MoleculeSystem *ms, CalcParams *params, double Temperature, double *m, double *q);
#endif // USE_MPI
// ----------------------------------------------------------
 
    
int assert_float_is_equal_to(double estimate, double true_value, double abs_tolerance);
 
double* linspace(double start, double end, size_t n);
double* arena_linspace(Arena *a, double start, double end, size_t n); 
size_t* arena_linspace_size_t(Arena *a, size_t start, size_t end, size_t n);
 
// ----------------------------------------------------------
// Histogram manipulation
// ----------------------------------------------------------
void send_histogram_and_reset(gsl_histogram *h);
void recv_histogram_and_append(Arena *a, int source, gsl_histogram **h);
gsl_histogram* gsl_histogram_extend_right(gsl_histogram* h, size_t add_bins);
int write_histogram_ext(FILE *fp, gsl_histogram *h, int count);
// ----------------------------------------------------------
 
// ----------------------------------------------------------
// String manipulation
// ----------------------------------------------------------
/* 
 * This struct represents a resizable buffer designed to build strings 
 * dynamically. It stores characters in a contiguous block of memory, 
 * allowing for manipulation of strings.
 */
typedef struct {
    char *items;
    size_t count;
    size_t capacity;
} String_Builder;
 
void sb_reserve(String_Builder *sb, size_t n);
void sb_append(String_Builder *sb, const char *line, size_t n);
void sb_append_null(String_Builder *sb);
void sb_append_cstring(String_Builder *sb, const char *line);
void sb_append_format(String_Builder *sb, const char *format, ...);
void sb_append_seconds_as_datetime_string(String_Builder *sb, int seconds);
void sb_reset(String_Builder *sb);
void sb_free(String_Builder *sb);
// ----------------------------------------------------------
 
// ----------------------------------------------------------
// Processing the results   
// ----------------------------------------------------------
CFnc copy_cfnc(CFnc cf);
SFnc copy_sfnc(SFnc sf);
Spectrum copy_spectrum(Spectrum sp); 
 
void normalize_cfnc(CFnc *cf);
 
void free_cfnc(CFnc cf);
void free_cfnc_array(CFncArray ca);
void free_sfnc(SFnc cf);
void free_spectrum(Spectrum sp); 
 
bool writetxt(const char *filename, double *x, double *y, size_t len, const char *header); 
 
bool write_correlation_function(const char *filename, CFnc cf);
int write_correlation_function_ext(FILE *fp, CFnc cf);
 
bool write_spectral_function(const char *filename, SFnc sf);
int write_spectral_function_ext(FILE *fp, SFnc sf);
 
bool write_spectrum(const char *filename, Spectrum sp);
int write_spectrum_ext(FILE *fp, Spectrum sp);
 
bool parse_number_of_trajectories_from_header(const char *header, double *ntraj);
bool parse_temperature_from_header(const char *header, double *Temperature); 
 
bool read_correlation_function(const char *filename, String_Builder *sb, CFnc *cf); 
bool read_spectral_function(const char *filename, String_Builder *sb, SFnc *sf); 
bool read_spectrum(const char *filename, String_Builder *sb, Spectrum *sp);
 
bool average_correlation_functions(CFnc *average, CFncs cfncs);
#define average_correlation_functions_ext(average, ...) average_correlation_functions__impl(average, sizeof((CFnc[]){__VA_ARGS__}) / sizeof(CFnc), __VA_ARGS__)
int average_correlation_functions__impl(CFnc *average, int arg_count,  ...);
 
void connes_apodization(Array a, double sampling_time); 
double *dct(double *v, size_t len);
double *idct(double *v, size_t len);
CFnc dct_sf_to_cf(SFnc sf);
SFnc idct_cf_to_sf(CFnc cf);
 
Spectrum compute_alpha(SFnc sf);
 
SFnc desymmetrize_d1(SFnc sf); 
SFnc desymmetrize_d2_sf(SFnc sf); 
Spectrum desymmetrize_d2_sp(Spectrum sp); 
SFnc desymmetrize_schofield_sf(SFnc sf); 
Spectrum desymmetrize_schofield_sp(Spectrum sp); 
SFnc desymmetrize_egelstaff(SFnc sf);
SFnc desymmetrize_egelstaff_from_cf(CFnc cf);
SFnc desymmetrize_frommhold(SFnc sf);
SFnc desymmetrize_frommhold_from_cf(CFnc cf);
CFnc egelstaff_time_transform(CFnc cf, bool frommhold_renormalization); 
Spectrum inv_desymmetrize_schofield(Spectrum sp);
Spectrum inv_desymmetrize_d2(Spectrum sp);
 
double* pad_to_power_of_two(double* v, size_t len, size_t* padded_len); 
 
extern WingParams INIT_WP;
double wingmodel(WingParams *wp, double t);
double wingmodel_image(WingParams *wp, double nu);
 
int wingmodel_f(const gsl_vector* x, void* data, gsl_vector* f);
int wingmodel_df(const gsl_vector* x, void* data, gsl_matrix * J);
 
void gsl_nonlinear_opt(size_t n, double* x, double* y, WingParams *wing_params);
WingParams fit_baseline(CFnc *cf, size_t EXT_RANGE_MIN);
 
 
/* Uses bit hack taken from: http://www.graphics.stanford.edu/~seander/bithacks.html#DetermineIfPowerOf2 */ 
static inline bool is_power_of_two(size_t n) {
    return n && !(n & (n - 1));
}
 
static inline unsigned int round_to_next_power_of_two(unsigned int n)
/*
 * rounds up to the next highest power of 2 using a clever bit hack
 * Taken from: 
 *  https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
 * 
 * works on 32-bit numbers
 */ 
{
    n--;
    n |= n >> 1;
    n |= n >> 2;
    n |= n >> 4;
    n |= n >> 8;
    n |= n >> 16;
    n++;
 
    return n;
}
    
#ifdef __cplusplus
}
#endif
 
#endif // HAWAII_H_