World.h

Go to the documentation of this file.

// 
//  This file is part of the FFEA simulation package
//  
//  Copyright (c) by the Theory and Development FFEA teams,
//  as they appear in the README.md file. 
// 
//  FFEA 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 3 of the License, or
//  (at your option) any later version.
// 
//  FFEA 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.
// 
//  You should have received a copy of the GNU General Public License
//  along with FFEA.  If not, see <http://www.gnu.org/licenses/>.
// 
//  To help us fund FFEA development, we humbly ask that you cite 
//  the research papers on the package.
//

#ifndef WORLD_H_INCLUDED
#define WORLD_H_INCLUDED

#include <stdio.h>
#include <stdlib.h>
#include <cstring>
#include <string>
#include <time.h>
#include <unistd.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <omp.h>
#include <ctime>
#include <algorithm>

#include <boost/algorithm/string.hpp>
#include <typeinfo>
#include <Eigen/Sparse>
#include <Eigen/Eigenvalues>

#ifdef FFEA_PARALLEL_FUTURE
#include <future>
#include <chrono>
#endif

// #include "MersenneTwister.h"
#include "RngStream.h"
#include "NearestNeighbourLinkedListCube.h"
#include "BEM_Poisson_Boltzmann.h"
#include "BiCGSTAB_solver.h"
#include "FFEA_user_info.h"
#include "FFEA_return_codes.h"
#include "FFEA_input_reader.h"
#include "mat_vec_types.h"
#include "mesh_node.h"
#include "tetra_element_linear.h"
#include "SimulationParams.h"
#include "Solver.h"
#include "SparseSubstitutionSolver.h"
#include "Face.h"
#include "Blob.h"
#include "World.h"
#include "VdW_solver.h"
#include "Steric_solver.h"
#include "LJSteric_solver.h"
#include "GenSoftSSINT_solver.h"
#include "PreComp_solver.h"
#include "LJ_matrix.h"
#include "BindingSite.h"
#include "Spring.h"
#include "SparseMatrixFixedPattern.h"
#include "KineticState.h"
#include "rod_structure.h"
#include "rod_blob_interface.h"

#include "dimensions.h"
using namespace std;

class World {
friend struct ffea_test; //allow  our unit test class to see ffea_world's private 
public:
    World();

    ~World();

    /* */
    int init(string FFEA_script_filename, int frames_to_delete, int mode, bool writeEnergy);

    /* */
    int get_smallest_time_constants();

    /* */
    int lem(set<int> blob_indices, int num_modes);

    /* */
    int dmm(set<int> blob_indices, int num_modes);

    /* */
    int dmm_rp(set<int> blob_indices, int num_modes);

    /* */
    int run();

    /* */
    int read_and_build_system(vector<string> script_vector);

    /* */
    int load_kinetic_maps(vector<string> map_fnames, vector<int> map_from, vector<int> map_to, int blob_index);

    /* */
    int build_kinetic_identity_maps();

    /* */
    int load_kinetic_states(string states_fname, int blob_index);

    /* */
    int load_kinetic_rates(string rates_fname, int blob_index);

    /* */
    void print_kinetic_rates_to_screen(int type);

    /* */
    void get_system_CoM(vector3 *system_CoM);

    /* */
    void get_system_centroid(vector3 *centroid);

    /* */
    void get_system_dimensions(vector3 *dimenstion_vector);

    /* */
    int enm(int *blob_index, int num_modes);

    /* */
    int get_num_blobs();


private:

    Blob **blob_array;

    Blob **active_blob_array;
    
    rod::Rod **rod_array;
    
    rod::Rod_blob_interface **rod_blob_interface_array;

    SparseMatrixFixedPattern ***kinetic_map;
    SparseMatrixFixedPattern ****kinetic_return_map;


    KineticState **kinetic_state;
    scalar ***kinetic_rate;
    scalar ***kinetic_base_rate;

    Spring *spring_array;

    int num_springs;

    int num_binding_sites;

    bool mass_in_system;

    int num_threads;

    RngStream *rng;

    unsigned long **Seeds;

    int num_seeds;

    RngStream *kinetic_rng;

    SimulationParams params;

    NearestNeighbourLinkedListCube lookup;

    FILE *trajectory_out;

    FILE *kinetics_out;

    FILE *measurement_out;

    bool writeDetailed;
    FILE *detailed_meas_out;

    FILE *trajbeads_out; 

    FFEA_input_reader *ffeareader;
    FFEA_input_reader *systemreader;
    

    scalar kineticenergy, strainenergy, springenergy, **springfieldenergy, ssintenergy, preCompenergy;

    vector3 L;


    vector3 CoM, CoG;
    scalar rmsd;

    FILE *checkpoint_out;

    /*
     *
     */
    //      SurfaceElementLookup surface_element_lookup;

    BEM_Poisson_Boltzmann PB_solver;

    int total_num_surface_faces;

    scalar *phi_Gamma;
    scalar *J_Gamma;
    scalar *work_vec;

    BiCGSTAB_solver nonsymmetric_solver;

    VdW_solver *vdw_solver;

    SSINT_matrix ssint_matrix;

    BindingSite_matrix binding_matrix;


    PreComp_params pc_params;
    PreComp_solver pc_solver;


    vector3 box_dim;

    long long step_initial;

    int load_springs(const char *fname);
    
    rod::Rod_blob_interface* rod_blob_interface_from_block(vector<string> block, int interface_id, FFEA_input_reader* systemreader, rod::Rod** rod_array, Blob** blob_array);

    rod::Rod* rod_from_block(vector<string> block, int block_id, FFEA_input_reader* systemreader);

    void activate_springs();

    int apply_springs();

    scalar get_spring_field_energy(int index0, int index1);

    int calculate_kinetic_rates();

    int choose_new_kinetic_state(int blob_index, int *target);

    int change_kinetic_state(int blob_index, int target_state);

    int get_next_script_tag(FILE *in, char *buf);

    void apply_dense_matrix(scalar *y, scalar *M, scalar *x, int N);

    void do_es();

    void make_trajectory_from_eigenvector(string traj_out_fname, int blob_index, int mode_index, Eigen_VectorX evec, scalar step);

    void print_evecs_to_file(string fname, Eigen_MatrixX ev, int num_rows, int num_modes);

    void print_evals_to_file(string fname, Eigen_VectorX ev, int num_modes, scalar scale);

    void write_eig_to_files(scalar *evals_ordered, scalar **evecs_ordered, int num_modes, int num_nodes);

    void write_output_header(FILE *fout, string fname);

    void print_trajectory_and_measurement_files(int step, scalar wtime);
    void print_checkpoints();
    void write_pre_print_to_trajfile(int step);
    void do_nothing(); 

    int prebuild_nearest_neighbour_lookup_wrapper(scalar cell_size); 
#ifdef FFEA_PARALLEL_FUTURE
    std::future<void> thread_writingTraj; 
    std::future<int> thread_updatingVdWLL; 
    std::future<int> thread_updatingPCLL; 
    bool updatingVdWLL(); 
    bool updatingVdWLL_ready_to_swap(); 
    int catch_thread_updatingVdWLL(int step, scalar wtime, int where); 
    bool updatingPCLL(); 
    bool updatingPCLL_ready_to_swap(); 
    int catch_thread_updatingPCLL(int step, scalar wtime, int where); 
#endif

    void make_measurements();

    void write_measurements_to_file(FILE *fout, int step);

    void write_detailed_measurements_to_file(FILE *fout);

    void print_trajectory_conformation_changes(FILE *fout, int step, int *from_index, int *to_index);

    void print_kinetic_files(int step);

    void print_static_trajectory(int step, scalar wtime, int blob_index);

    void calc_blob_corr_matrix(int num_blobs,scalar *blob_corr);

    scalar *blob_corr;

    int die_with_dignity(int step, scalar wtime); 
};

#endif