ForwardFlux.h

 
#pragma once 
 
#include "Method.h"
#include <fstream>
#include <random>
#include <deque>
#include "sys/stat.h"
 
namespace SSAGES
{
 
    class ForwardFlux : public Method
    {
    protected:
        class FFSConfigID
        {
           public:
            unsigned int l;     
            unsigned int n;     
            unsigned int a;     
            unsigned int lprev; 
            unsigned int nprev; 
            unsigned int aprev; 
 
 
            FFSConfigID(const unsigned int l, 
                        const unsigned int n, 
                        const unsigned int a, 
                        const unsigned int lprev, 
                        const unsigned int nprev, 
                        const unsigned int aprev): 
             l(l),n(n),a(a),lprev(lprev),nprev(nprev),aprev(aprev)
            {}
 
            FFSConfigID():
             l(0),n(0),a(0),lprev(0),nprev(0),aprev(0)
            {}
        };
 
        double _ninterfaces;
 
        std::vector<double> _interfaces;
 
        bool _interfaces_increase;
 
        double _cvvalue_previous;
 
        double _cvvalue;
 
        double _rate;
 
        std::vector<FFSConfigID> Lambda0ConfigLibrary;
 
        double _N0TotalSimTime;
        
        unsigned int _N0Target;
 
        double _fluxA0;
 
        std::vector<unsigned int> _M;
 
        std::vector<unsigned int> _A;
 
        //bool _computefluxA0;
 
        std::vector<double> _P;
 
        std::vector<unsigned int> _S;
 
        std::vector<unsigned int> _N ;
        
        bool _pop_tried_but_empty_queue;
 
        bool _initialFluxFlag;
 
        bool initializeQueueFlag;
 
        FFSConfigID myFFSConfigID;
 
        bool _saveTrajectories;
 
        unsigned int _nfailure_total;
 
        std::vector<std::vector<double>> _pB;
 
        unsigned int _current_interface;
 
        std::deque<FFSConfigID> FFSConfigIDQueue; 
 
        std::string _output_directory;
 
        std::ofstream _trajectory_file;
        
        std::default_random_engine _generator;
 
        unsigned int iteration_;
 
        //-----------------------------------------------------------------
        // Private Functions
        //-----------------------------------------------------------------
        
        void CheckInitialStructure(const CVList&);
 
        void WriteInitialFlux();
 
        void AddNewIDsToQueue();
 
 
        bool HasReturnedToA(double current);
 
 
        int HasCrossedInterface(double current, double prev, unsigned int i);
 
 
        void WriteFFSConfiguration(Snapshot *snapshot, FFSConfigID& ffsconfig, bool wassuccess);
 
 
        void ReadFFSConfiguration(Snapshot *snapshot, FFSConfigID& ffsconfig, bool wassuccess);
 
        void ComputeInitialFlux(Snapshot*, const CVList&);
 
        virtual void CheckForInterfaceCrossings(Snapshot*, const class CVManager&) = 0;
 
        virtual void InitializeQueue(Snapshot*, const CVList&) =0;
 
 
        void ComputeTransitionProbabilities();
 
        void PrintQueue();
 
        void PopQueueMPI(Snapshot*, const CVList&, unsigned);
        
 
        void FluxBruteForce(Snapshot*, const CVList&);
        
        void ReconstructTrajectories(Snapshot *);
 
        void ComputeCommittorProbability(Snapshot *);
        
        void AppendTrajectoryFile(Snapshot*, std::ofstream&);
 
        void OpenTrajectoryFile(std::ofstream&);
 
    public:
 
        ForwardFlux(const MPI_Comm& world,
                    const MPI_Comm& comm, 
                    double ninterfaces, std::vector<double> interfaces,
                    unsigned int N0Target, std::vector<unsigned int> M,
                    bool initialFluxFlag, bool saveTrajectories,
                    unsigned int currentInterface, std::string output_directory, unsigned int frequency) : 
         Method(frequency, world, comm), _ninterfaces(ninterfaces), _interfaces(interfaces), _N0Target(N0Target), 
         _M(M), _initialFluxFlag(initialFluxFlag), _saveTrajectories(saveTrajectories), _current_interface(currentInterface),
         _output_directory(output_directory), _generator(1), iteration_(0) 
        {
            //_output_directory = "FFSoutput";
            mkdir(_output_directory.c_str(),S_IRWXU); //how to make directory?
 
            //_current_interface = 0;
            _N0TotalSimTime = 0;
 
            if (!_initialFluxFlag)
                initializeQueueFlag = true;
 
            _A.resize(_ninterfaces);
            _P.resize(_ninterfaces);
            _S.resize(_ninterfaces);
            _N.resize(_ninterfaces);
 
            //_N[_current_interface] = _N0Target;
            /*_M.resize(_ninterfaces);
 
            //code for setting up simple simulation and debugging
            _ninterfaces = 5;
            int i;
            _interfaces.resize(_ninterfaces);
            _interfaces[0]=-1.0;
            _interfaces[1]=-0.95;
            _interfaces[2]=-0.8;
            _interfaces[3]= 0;
            _interfaces[4]= 1.0;
 
            _saveTrajectories = true;
 
            _initialFluxFlag = true;
 
            for(i=0;i<_ninterfaces;i++) _M[i] = 50;
 
            _N0Target = 100;*/
            //_N0Target = 100;
            _nfailure_total = 0;
 
 
            //check if interfaces monotonically increase or decrease
            bool errflag = false;
            if (_interfaces[0]< _interfaces[1]) _interfaces_increase = true;
            else if (_interfaces[0] > _interfaces[1]) _interfaces_increase = false;
            else errflag = true;
            for (unsigned int i=0;i<_ninterfaces-1;i++)
            {
                if((_interfaces_increase) && (_interfaces[i] >= _interfaces[i+1]))
                    errflag = true;
                else if ((!_interfaces_increase) && (_interfaces[i] <= _interfaces[i+1]))
                    errflag = true;
            }
            if (errflag)
            {
                std::cerr << "Error! The interfaces are poorly defined. They must be monotonically increasing or decreasing and cannot equal one another! Please fix this.\n";
                for (auto interface : _interfaces)
                    std::cerr << interface << " ";
                std::cerr << "\n";
                MPI_Abort(world_, EXIT_FAILURE);
            }
 
            // This is to generate an artificial Lambda0ConfigLibrary, Hadi's code does this for real
            // THIS SHOULD BE SOMEWHERE ELSE!!! 
            Lambda0ConfigLibrary.resize(_N0Target);
            std::normal_distribution<double> distribution(0,1);
            for (unsigned int i = 0; i < _N0Target ; i++)
            {
                Lambda0ConfigLibrary[i].l = 0;
                Lambda0ConfigLibrary[i].n = i;
                Lambda0ConfigLibrary[i].a = 0;
                Lambda0ConfigLibrary[i].lprev = 0;
                Lambda0ConfigLibrary[i].nprev = i;
                Lambda0ConfigLibrary[i].aprev = 0;
                //FFSConfigID ffsconfig = Lambda0ConfigLibrary[i];
            }
            /*// Write the dump file out
            std::ofstream file;
            std::string filename = _output_directory + "/l" + std::to_string(ffsconfig.l) + "-n" + std::to_string(ffsconfig.n) + ".dat";
            file.open(filename.c_str());
 
            //first line gives ID of where it came from
            file << ffsconfig.lprev << " " << ffsconfig.nprev << " " << ffsconfig.aprev << "\n";
            //write position and velocity
            file << "1 -1 0 0 " << distribution(_generator) << " "  << distribution(_generator) << " 0\n";
 
            }
            */
            _pop_tried_but_empty_queue = false;
        }
 
        void PreSimulation(Snapshot* snapshot, const class CVManager& cvmanager) override;
 
        virtual void PostIntegration(Snapshot* snapshot, const class CVManager& cvmanager) override = 0;
 
        void PostSimulation(Snapshot* snapshot, const class CVManager& cvmanager) override;
 
        static ForwardFlux* Build(const Json::Value& json, 
                                  const MPI_Comm& world,
                                  const MPI_Comm& comm,
                                  const std::string& path);
    };
}
 
/*
File Formats:
_indexfile
interface(some integer) dump_file_name(a string that contains interface and trial number)
example: 1 dump_1_10.xyz
 
dumpfile
atomid posx posy posz vx vy vz
 
 
*/