Atomistic simulation classes

atomistic.C

// Class forward declarations
class atomisticMoleculeCloud;

/*---------------------------------------------------------------------------*\
                           Class atomisticMolecule Declaration
\*---------------------------------------------------------------------------*/

class atomisticMolecule
:
    public Particle<atomisticMolecule>
{

public:

    // Values of special that are less than zero are for built-in functionality.
    // Values greater than zero are user specifiable/expandable (i.e. test
    // special_ >= SPECIAL_USER)

    enum specialTypes
    {
        SPECIAL_TETHERED = -1,
        SPECIAL_FROZEN   = -2,
        NOT_SPECIAL      = 0,
        SPECIAL_USER     = 1
    };

    //- Class to hold molecule constant properties
    class constantProperties
    {
        // Private data

            constPropSite site_;

            scalar mass_;

        // Private Member Functions

    public:

        inline constantProperties();

        //- Construct from dictionary
        inline constantProperties
        (
            const dictionary& dict,
            const reducedUnits& rU,
            const label& siteId
        );

        // Member functions

        inline const constPropSite& site() const;

        inline label degreesOfFreedom() const;

//         inline const bool& pairPotentialSite() const;
// 
//         inline const bool& electrostaticSite() const;

        inline scalar mass() const;

    };


    //- Class used to pass tracking data to the trackToFace function
    class trackData
    :
        public Particle<atomisticMolecule>::trackData
    {
        atomisticMoleculeCloud& molCloud_;

        // label specifying which part of the integration algorithm is taking
        label part_;


    public:

        // Constructors

            trackData
            (
                atomisticMoleculeCloud& molCloud,
                label part
            );

        // Member functions

            inline atomisticMoleculeCloud& molCloud();

            inline label part() const;
    };


private:

    // Private data

        //- Be careful with the ordering of data.
        //  It has an impact on binary transfer:
        //    -# Put the largest data members 1st
        //    -# Pair up labels,
        //    -# Don't go scalar-label, scalar-label, because in 64bit mode,
        //       the labels will be padded by 4bytes.

        //- linear velocity
        vector v_;

        //- linear acceleration
        vector a_;

        //- force
        vector f_;

        vector specialPosition_;

        scalar potentialEnergy_;

        // - r_ij outer product f_ij: virial contribution
//         tensor rDotf_;
        tensor rf_;

        label special_;

        label id_;

        //- distance to the closest molecule
        scalar R_;

        scalar frac_;

        //- each molecule will have a unique label
        label trackingNumber_;


        //- is molecule in a collision
        bool collision_;


public:

    friend class Cloud<atomisticMolecule>;

    // Constructors

        //- Construct from components
        inline atomisticMolecule
        (
            const Cloud<atomisticMolecule>& c,
            const vector& position,
            const label celli,
            const vector& v,
            const vector& f,
            const vector& specialPosition,
            const constantProperties& constProps,
            const scalar& pE,
            const tensor& rf,
            const label special,
            const label id,
            const label trackingNumber,
            const scalar& frac,
            const bool collision
        );

        //- Construct from Istream
        atomisticMolecule
        (
            const Cloud<atomisticMolecule>& c,
            Istream& is,
            bool readFields = true
        );

        //- Construct and return a clone
        autoPtr<atomisticMolecule> clone() const
        {
            return autoPtr<atomisticMolecule>(new atomisticMolecule(*this));
        }


    // Member Functions

        // Tracking

            bool move(trackData&);

            void transformProperties(const tensor& T);

            void transformProperties(const vector& separation);

//             void setSiteSizes(label size);


        // Access


            inline const vector& v() const;
            inline vector& v();

            inline const vector& a() const;
            inline vector& a();

            inline const vector& f() const;
            inline vector& f();
//             inline const List<vector>& siteForces() const;
//             inline List<vector>& siteForces();
// 
//             inline const List<vector>& sitePositions() const;
//             inline List<vector>& sitePositions();

            inline const vector& specialPosition() const;
            inline vector& specialPosition();

            inline scalar potentialEnergy() const;
            inline scalar& potentialEnergy();

            // - Return virial contribution
            inline const tensor& rf() const;
            inline tensor& rf();

            inline label special() const;
            inline label& special();

            inline bool tethered() const;

            inline label id() const;
            inline label& id();


            // - return shortest radius
            inline const scalar& R() const;
            inline scalar& R();

            // - return tracking number
            inline const label& trackingNumber() const;
            inline label& trackingNumber();

            // - return fraction number
            inline const scalar& fraction() const;
            inline scalar& fraction();

            // - return collision status
            inline const bool& collisionStatus() const;
            inline bool& collisionStatus();


//             inline const bool& move() const;
//             inline bool& move();

            inline bool frozen();

    // Member Operators

        //- Overridable function to handle the particle hitting a patch
        //  Executed before other patch-hitting functions
        bool hitPatch
        (
            const polyPatch&,
            atomisticMolecule::trackData& td,
            const label patchI
        );

        //- Overridable function to handle the particle hitting a patch
        //  Executed before other patch-hitting functions without trackData
        bool hitPatch
        (
            const polyPatch& p,
            int& td,
            const label patchI
        );

        //- Overridable function to handle the particle hitting a processorPatch
        void hitProcessorPatch
        (
            const processorPolyPatch&,
            atomisticMolecule::trackData& td,
            const label& patchi
        );

        //- Overridable function to handle the particle hitting a processorPatch
        //  without trackData
        void hitProcessorPatch
        (
            const processorPolyPatch&,
            int&,
            const label& patchi
        );

        //- Overridable function to handle the particle hitting a wallPatch
        void hitWallPatch
        (
            const wallPolyPatch&,
            atomisticMolecule::trackData& td
        );

        //- Overridable function to handle the particle hitting a wallPatch
        //  without trackData
        void hitWallPatch
        (
            const wallPolyPatch&,
            int&
        );

        //- Overridable function to handle the particle hitting a polyPatch
        void hitPatch
        (
            const polyPatch&,
            atomisticMolecule::trackData& td
        );

        //- Overridable function to handle the particle hitting a polyPatch
        //  without trackData
        void hitPatch
        (
            const polyPatch&,
            int&
        );

        void hitCyclicPatch
        (
            const cyclicPolyPatch& cpp,
            int&,
            const label& patchi
        );

        void hitCyclicPatch
        (
            const cyclicPolyPatch& cpp,
            atomisticMolecule::trackData& td,
            const label& patchi
        );

        static void readFields(Cloud<atomisticMolecule>& mC);

        static void writeFields(const Cloud<atomisticMolecule>& mC);


    // IOstream Operators

        friend Ostream& operator<<(Ostream&, const atomisticMolecule&);
};

class atomisticMoleculeCloud
:
    public Cloud<atomisticMolecule>
{

private:

    // Private data

        const polyMesh& mesh_;

        const potential& pot_;

        const reducedUnits& redUnits_;

        List<DynamicList<atomisticMolecule*> > cellOccupancy_;

//         atomisticInteractionLists il_;

        List<atomisticMolecule::constantProperties> constPropList_;

        //- measurement fields
        atomisticFieldProperties fields_; 

        // cyclic boundaries
        cyclicBoundaries cyclics_;

        //- boundaries
        atomisticBoundaries boundariesInfo_;

        //- controllers
        atomisticControllers controllers_;

        //- atomisticMolecule tracking
        atomisticFaceTracker trackingInfo_;

        //- index which is used for assigining tracking numbers
        //  to newly created molecules.
        trackingNumber moleculeTracking_;

//         atomisticKernel interactions_;

        mdKernel<atomisticMolecule> mdKernel_;

        List<DynamicList<atomisticMolecule*> > interactionList_;

//         clockTimer forceComputation_;
//         clockTimer forcePreparation_;
//         clockTimer move_;
        
        // required for mean free path measurements
        DynamicList<label> collidedMolecules_;
        DynamicList<vector> collidedMoleculesPosition_; 
        
        Random rndGen_;


    // Private Member Functions

        void buildConstProps();

//         void setSiteSizesAndPositions();

        //- Determine which molecules are in which cells

//         void buildCellOccupancyNoInteractionLists();

//         void calculatePairForce();

//         void calculateTetherForce();

        void removeHighEnergyOverlaps();

        inline vector equipartitionLinearVelocity
        (
            scalar temperature,
            scalar mass
        );



        void setTrackingNumbers();

        void computePairForces();
        void computePairForces2();
        void accelerationUpdate();
        void velocityUpdate(const scalar& trackTime);

        //- Disallow default bitwise copy construct
        atomisticMoleculeCloud(const atomisticMoleculeCloud&);

        //- Disallow default bitwise assignment
        void operator=(const atomisticMoleculeCloud&);


public:

    // Static data members

    // Constructors

        //- Construct for mdFoam
        
        atomisticMoleculeCloud
        (
            Time& t,
            const polyMesh& mesh,
            const potential& pot,
            const reducedUnits& rU
        );
        
        atomisticMoleculeCloud
        (
            Time& t,
            const word& cloudName,         
            const polyMesh& mesh,
            const potential& pot,
            const reducedUnits& rU
        );


        //- General constructor
        atomisticMoleculeCloud
        (
            Time& t,
            const polyMesh& mesh,
            const potential& pot,
            const reducedUnits& rU,
            const word& option,
            const bool& clearFields
        );


    // Static constructors

        static autoPtr<atomisticMoleculeCloud> New
        (
            Time& t,
            const polyMesh& mesh,
            const potential& pot,
            const reducedUnits& rU
        );

        static autoPtr<atomisticMoleculeCloud> New
        (
            Time& t,
            const word& cloudName,          
            const polyMesh& mesh,
            const potential& pot,
            const reducedUnits& rU
        );
        
        static autoPtr<atomisticMoleculeCloud> New
        (
            Time& t,
            const polyMesh& mesh,
            const potential& pot,
            const reducedUnits& rU,
            const word& option,
            const bool& clearFields
        );


    // Member Functions

        //- Evolve the molecules (move, calculate forces, control state etc)
        void preliminaries();
        void initialHalfVelocity();
        void move();
        void controlBeforeForces();
        void updateAcceleration();
        void controlAfterForces();
        void finalHalfVelocity();
        void postPreliminaries();


        void evolve();
        void evolveBeforeForces();
        void evolveAfterForces();
//         void evolveBeforeForcesGPU();

        void clearLagrangianFields();

        void calculateForce();

        void setMoleculesForDirectInteraction();

        void checkMoleculesInMesh();
        void buildCellOccupancy();
        void rebuildCellOccupancy();

        void prepareKernel();

        // insert a molecule in the moleculeCloud
        void createMolecule
        (
            const vector& position,
            const label cell,
            const vector& v,
            const vector& f,
            const vector& specialPosition,
            const scalar& pE,
            const tensor& rDotf,
            const label special,
            const label id,
            const label trackingNumber,
            const scalar& fraction,
            const bool collision
        );
};