action.h

Go to the documentation of this file.

 
#include "constants.h"
 
#ifndef ACTION_H 
#define ACTION_H
 
class Path;
class PotentialBase;
class LookupTable;
class WaveFunctionBase;
 
// ========================================================================  
// ActionBase Class
// ========================================================================  
class ActionBase {
 
    public:
        ActionBase (const Path &, LookupTable &, PotentialBase *, 
                PotentialBase *, WaveFunctionBase *, bool _local=true,
                string _name="Base", double _endFactor=1.0, int _period=1);
        virtual ~ActionBase();
 
        string getActionName () { return name; }
 
        double kineticAction ();
        double kineticAction (const beadLocator &);
        double kineticAction (const beadLocator &, int wlLength);
 
        virtual double potentialAction() {return 0.0;}
        virtual double potentialAction (const beadLocator &, const beadLocator &);
        virtual double potentialAction (const beadLocator &) { return 0.0; }
 
        /* The bare potential action and its correction */
        virtual double barePotentialAction (const beadLocator &) { return 0.0; }
        virtual double potentialActionCorrection (const beadLocator &) {return 0.0; }
        virtual double potentialActionCorrection (const beadLocator &, const beadLocator &) { return 0.0; }
 
        /* Various derivatives of the potential action */
        virtual double derivPotentialActionTau (int) { return 0.0; }
        virtual double derivPotentialActionLambda (int) { return 0.0; }
        virtual double secondderivPotentialActionTau (int) { return 0.0; }
 
        /* Various derivatives of the potential action with a cutoff */
        virtual double derivPotentialActionTau (int,double) { return 0.0; }
        virtual double derivPotentialActionLambda (int,double) { return 0.0; }
        
        /* gradient of the potential action */
        virtual dVec gradPotentialAction(int) { return 0.0; }
 
        /* r \dot gradU -> for virial estimator.
         * It is necessary to split the terms in order to compute
         * the specific heat efficiently.*/
        virtual double rDOTgradUterm1(int) { return 0.0; }
        virtual double rDOTgradUterm2(int) { return 0.0; }
 
        /* (r-R) \dot gradU -> for centroid virial estimator.
         * It is necessary to split the terms in order to compute
         * the specific heat efficiently.*/
        virtual double deltaDOTgradUterm1(int) { return 0.0; }
        virtual double deltaDOTgradUterm2(int) { return 0.0; }
 
        /* return virial kinetic energy estimator term */
        virtual double virKinCorr(int) { return 0.0; }
 
        /* The bare local potential at a single time slice */
        virtual TinyVector<double,2> potential(int) { return TinyVector<double,2>(0.0); }
        virtual double potential(int,double) { return 0.0; }
 
        void setShift(int _shift) { shift = _shift; }
        
        int getShift() { return shift; }
 
        double rho0(const dVec&, const dVec&, int); 
        double rho0(const beadLocator&, const beadLocator&, int); 
        double rho0(const dVec&, const int); 
        
        double ensembleWeight(const int);
 
        const bool local;
 
        /* The period of the action */
        const int period;
 
        PotentialBase *externalPtr;     
        PotentialBase *interactionPtr;  
 
        Array <int,1> sepHist;          
        Array <int,1> cylSepHist;       
 
    protected:
        string name;                    
 
        LookupTable &lookup;            
        const Path &path;               
        WaveFunctionBase *waveFunctionPtr;   
        double endFactor;               
 
        int shift;                      
 
        /* These definitions are needed for weighting in the canonical ensemble */
        bool canonical;                 
        int numBeads0;                  
        bool window;                    // Whether or not to force particle numberin to window
        int windowWidth;                // Half width of particle number window
        bool gaussianEnsemble;          // Whether or not to use gaussian ensemble weight
        double gaussianEnsembleSD;        // Standard deviation of guassian ensemble weight
 
        double tau() {return shift * constants()->tau();}
 
        beadLocator bead2,bead3;        // Bead indexers
        dVec sep,sep2;                  // The spatial separation between beads.
        double dSep;                    // The discretization for the separation histogram
 
        /* Update the separation histogram */
        void updateSepHist(const dVec &);   
};
 
// ========================================================================  
// LocalAction Class
// ========================================================================  
class LocalAction : public ActionBase {
 
    public:
        LocalAction (const Path &, LookupTable &, PotentialBase *, 
                PotentialBase *, WaveFunctionBase *, const TinyVector<double,2>&, 
                const TinyVector<double,2>&, bool _local=true, string _name="Local",
                double _endFactor=1.0, int _period=1);
        virtual ~LocalAction();
 
        /* The potential action */
        double potentialAction ();
        double potentialAction (const beadLocator &);
 
        /* The bare potential action and its correction */
        double barePotentialAction (const beadLocator &);
        double potentialActionCorrection (const beadLocator &);
        double potentialActionCorrection (const beadLocator &, const beadLocator &);
 
        /* Various derivatives of the potential action */
        double derivPotentialActionTau (int);
        double derivPotentialActionLambda (int);
        double secondderivPotentialActionTau (int);
        double derivPotentialActionTau (int,double);
        double derivPotentialActionLambda (int,double);
 
        /* gradient of the potential action */
        virtual dVec gradPotentialAction(int slice) { return gradU(slice); }
 
        /* R \dot gradU -> for virial estimators */
        double rDOTgradUterm1(int); 
        double rDOTgradUterm2(int); 
 
        /* delta \dot gradU -> for centroid virial estimators */
        virtual double deltaDOTgradUterm1(int slice) { return deltadotgradUterm1(slice); }
        virtual double deltaDOTgradUterm2(int slice) { return deltadotgradUterm2(slice); }
 
        /* Returns virial kinetic energy correction term. */
        virtual double virKinCorr(int slice) { return virialKinCorrection(slice); }
 
        /* The bare local potential at a single time slice */
        virtual TinyVector<double,2> potential(int slice) { return V(slice); }
        virtual double potential(int slice, double maxR) { return V(slice,maxR); }
 
    protected:
        int eo;                         
 
        TinyVector <double,2> VFactor;      
        TinyVector <double,2> gradVFactor;  
 
        /* The full potential for a single bead and all beads at a single
         * time slice. */
        double V(const beadLocator&);   
        double V(const int, const double);
 
        /* For the potential at a given time slice we separate the interaction
         * and potential parts */
        TinyVector <double,2> V(const int);
 
        /* The gradient of the potential squared for a single bead and all beads
         * at a single time slice. */
        double gradVSquared(const beadLocator&);    
        double gradVSquared(const int);
        double gradVSquared(const int, const double);
 
        /* The full potential with the NN lookup table for a single bead and all
         * beads at a single time slice. */
        double Vnn(const beadLocator&);
        double Vnn(const int);
 
        /* The bare potential action for a trajectory */
        double bareUnn(const beadLocator &, const beadLocator &);
 
        /* The gradient of the potential squared for a single bead using the 
         * nearest neighbor lookup table */
        double gradVnnSquared(const beadLocator&);
 
        /* gradient and Laplacian of potential energy */
        dVec gradientV(const int);
 
        /* T-matrix needed for grad((gradV)^2) */
        dMat tMatrix(const int);
 
        /* r \dot gradU -> for virial estimator */
        double RdotgradUterm1(const int);
        double RdotgradUterm2(const int);
 
        /* delta \dot gradU -> for centroid virial estimator */
        double deltadotgradUterm1(const int);
        double deltadotgradUterm2(const int);
 
        /* virial kinetic energy term */
        double virialKinCorrection(const int);
 
        /* gradient of potential action */
        dVec gradU(const int);
 
};
 
// ========================================================================  
// NonLocalAction Class
// ========================================================================  
class NonLocalAction : public ActionBase {
 
    public:
        NonLocalAction (const Path &, LookupTable &, PotentialBase *, 
                PotentialBase *, WaveFunctionBase *, bool _local=false,
                string _name="Non-local");
        virtual ~NonLocalAction();
 
        /* The potential Action */
        double potentialAction ();
        double potentialAction (const beadLocator &);
 
        /* Derivatives of the potential action */
        double derivPotentialActionTau (int);
        double derivPotentialActionLambda (int);
    
        /* The bare local external and interaction potential at a single time slice */
        virtual TinyVector<double,2> potential(int slice) { return U(slice); }
 
    protected:
        TinyVector<double,2> U(int);
    
    private:
        vector<bool> NNbead; //Records which beads where already visited for NN operation
};
#endif