MdSystem.h

#ifndef MCMD_MD_SYSTEM_H
#define MCMD_MD_SYSTEM_H

/*
* Simpatico - Simulation Package for Polymeric and Molecular Liquids
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <mcMd/simulation/System.h>     // base class
#include <util/signal/Signal.h>         // members

#include <util/global.h>

namespace Util {
   template <typename T> class Factory;
}

namespace McMd
{

   using namespace Util;

   // Forward declarations

   #ifndef SIMP_NOPAIR
   class MdPairPotential;
   #endif
   #ifdef SIMP_BOND
   class BondPotential;
   #endif
   #ifdef SIMP_ANGLE
   class AnglePotential;
   #endif
   #ifdef SIMP_DIHEDRAL
   class DihedralPotential;
   #endif
   #ifdef SIMP_COULOMB
   class MdCoulombPotential;
   #endif
   #ifdef SIMP_EXTERNAL
   class ExternalPotential;
   #endif
   #ifdef SIMP_SPECIAL
   class SpecialPotential;
   #endif
   #ifdef SIMP_TETHER
   class TetherPotential;
   #endif

   class MdIntegrator;
   class McSystem;

   class MdSystem : public System
   {

   public:

      MdSystem();

      MdSystem(McSystem& system);

      virtual ~MdSystem();



      virtual void readParameters(std::istream &in);

      virtual void loadParameters(Serializable::IArchive &ar);

      virtual void saveParameters(Serializable::OArchive &ar);



      virtual void readConfig(std::istream& in);

      using System::readConfig;

      virtual void loadConfig(Serializable::IArchive& ar);

      void generateMolecules(Array<int> const & capacities,
                             Array<double> const & diameters);



      void calculateForces();

      double kineticEnergy() const;

      double potentialEnergy();

      void unsetPotentialEnergy();

      template <typename T>
      void computeStress(T& stress) const;

      template <typename T>
      void computeVirialStress(T& stress) const;

      void unsetVirialStress();

      template <typename T>
      void computeKineticStress(T& stress) const;



      #ifndef SIMP_NOPAIR

      MdPairPotential& pairPotential() const;
      #endif

      bool hasBondPotential() const;

      #ifdef SIMP_BOND

      BondPotential& bondPotential() const;
      #endif

      #ifdef SIMP_ANGLE

      bool hasAnglePotential() const;

      AnglePotential& anglePotential() const;
      #endif

      #ifdef SIMP_DIHEDRAL

      bool hasDihedralPotential() const;

      DihedralPotential& dihedralPotential() const;
      #endif

      #ifdef SIMP_COULOMB

      bool hasCoulombPotential() const;

      MdCoulombPotential& coulombPotential() const;
      #endif

      #ifdef SIMP_EXTERNAL

      bool hasExternalPotential() const;

      ExternalPotential& externalPotential() const;
      #endif

      #ifdef SIMP_SPECIAL

      bool hasSpecialPotential() const;

      SpecialPotential& specialPotential() const;
      #endif

      #ifdef MCMD_LINK

      bool hasLinkPotential() const;

      BondPotential& linkPotential() const;
      #endif

      #ifdef SIMP_TETHER

      TetherPotential& tetherPotential() const;
      #endif



      void setZeroForces();

      void setZeroVelocities();

      void setBoltzmannVelocities(double temperature);

      Vector removeDriftVelocity();

      void shiftAtoms();



      Signal<>& positionSignal();

      Signal<>& velocitySignal();



      MdIntegrator& mdIntegrator();

      Factory<MdIntegrator>& mdIntegratorFactory();

      virtual bool isValid() const;


   protected:

      virtual ConfigIo* newDefaultConfigIo();

   private:

      Signal<>  positionSignal_;

      Signal<>  velocitySignal_;

      #ifndef SIMP_NOPAIR
      MdPairPotential* pairPotentialPtr_;
      #endif

      #ifdef SIMP_BOND
      BondPotential* bondPotentialPtr_;
      #endif

      #ifdef SIMP_ANGLE
      AnglePotential* anglePotentialPtr_;
      #endif

      #ifdef SIMP_DIHEDRAL
      DihedralPotential* dihedralPotentialPtr_;
      #endif

      #ifdef SIMP_COULOMB
      MdCoulombPotential* coulombPotentialPtr_;
      #endif

      #ifdef SIMP_EXTERNAL
      ExternalPotential* externalPotentialPtr_;
      #endif

      #ifdef SIMP_SPECIAL
      SpecialPotential* specialPotentialPtr_;
      #endif

      #ifdef MCMD_LINK
      BondPotential* linkPotentialPtr_;
      #endif

      #ifdef SIMP_TETHER
      TetherPotential* tetherPotentialPtr_;
      #endif

      MdIntegrator *mdIntegratorPtr_;

      Factory<MdIntegrator> *mdIntegratorFactoryPtr_;

      bool createdMdIntegratorFactory_;

      /*
      * Implementations of the explicit specializations of the public
      * stress calculators computeVirialStress(T& ) etc. for T = double,
      * Util::Vector and Util::Tensor simply call these private method
      * templates. This allows a single private template method to be
      * used to compile the 3 relevant explicit specializations of each
      * stress calculator, while allowing the template implementations
      * to be defined in MdSystem.cpp, rather than in this header file.
      */

      template <typename T>
      void computeKineticStressImpl(T& stress) const;

      template <typename T>
      void computeVirialStressImpl(T& stress) const;

   };

   // Inline functions

   #ifndef SIMP_NOPAIR
   /*
   * Return PairPotential by reference.
   */
   inline MdPairPotential& MdSystem::pairPotential() const
   {
      assert(pairPotentialPtr_);
      return *pairPotentialPtr_;
   }
   #endif

   #ifdef SIMP_BOND
   /*
   * Does a bond potential exist?
   */
   inline bool MdSystem::hasBondPotential() const
   {  return bool(bondPotentialPtr_); }

   /*
   * Return bond potential by reference.
   */
   inline BondPotential& MdSystem::bondPotential() const
   {
      assert(bondPotentialPtr_);
      return *bondPotentialPtr_;
   }
   #endif

   #ifdef SIMP_ANGLE
   /*
   * Does an angle potential exist?
   */
   inline bool MdSystem::hasAnglePotential() const
   {  return bool(anglePotentialPtr_); }

   /*
   * Return angle potential by reference.
   */
   inline AnglePotential& MdSystem::anglePotential() const
   {
      assert(anglePotentialPtr_);
      return *anglePotentialPtr_;
   }
   #endif

   #ifdef SIMP_DIHEDRAL
   /*
   * Does a dihedral potential exist?
   */
   inline bool MdSystem::hasDihedralPotential() const
   {  return bool(dihedralPotentialPtr_); }

   /*
   * Return dihedral potential by reference.
   */
   inline DihedralPotential& MdSystem::dihedralPotential() const
   {
      assert(dihedralPotentialPtr_);
      return *dihedralPotentialPtr_;
   }
   #endif

   #ifdef SIMP_COULOMB
   /*
   * Does a Coulomb potential exist?
   */
   inline bool MdSystem::hasCoulombPotential() const
   {  return bool(coulombPotentialPtr_); }

   /*
   * Return Coulomb potential by reference.
   */
   inline MdCoulombPotential& MdSystem::coulombPotential() const
   {
      assert(coulombPotentialPtr_);
      return *coulombPotentialPtr_;
   }
   #endif

   #ifdef SIMP_EXTERNAL
   /*
   * Does an external potential exist?
   */
   inline bool MdSystem::hasExternalPotential() const
   {  return bool(externalPotentialPtr_); }

   /*
   * Return external potential by reference.
   */
   inline ExternalPotential& MdSystem::externalPotential() const
   {
      assert(externalPotentialPtr_);
      return *externalPotentialPtr_;
   }
   #endif

   #ifdef SIMP_SPECIAL
   /*
   * Does a special potential exist?
   */
   inline bool MdSystem::hasSpecialPotential() const
   {  return bool(specialPotentialPtr_); }

   /*
   * Return special potential by reference.
   */
   inline SpecialPotential& MdSystem::specialPotential() const
   {
      assert(specialPotentialPtr_);
      return *specialPotentialPtr_;
   }
   #endif

   #ifdef MCMD_LINK
   /*
   * Does a link potential exist?
   */
   inline bool MdSystem::hasLinkPotential() const
   {  return bool(linkPotentialPtr_); }

   /*
   * Return link potential by reference.
   */
   inline BondPotential& MdSystem::linkPotential() const
   {
      assert(linkPotentialPtr_);
      return *linkPotentialPtr_;
   }
   #endif

   #ifdef SIMP_TETHER
   /*
   * Return tether potential by reference.
   */
   inline TetherPotential& MdSystem::tetherPotential() const
   {
      assert(tetherPotentialPtr_);
      return *tetherPotentialPtr_;
   }
   #endif

   /*
   * Return the MdIntegrator by reference.
   */
   inline MdIntegrator& MdSystem::mdIntegrator()
   {
      assert(mdIntegratorPtr_);
      return *mdIntegratorPtr_;
   }

   /*
   * Signal to indicate change in atomic positions.
   */
   inline Signal<>& MdSystem::positionSignal()
   {  return positionSignal_; }

   /*
   * Signal to indicate change in atomic velocities.
   */
   inline Signal<>& MdSystem::velocitySignal()
   {  return velocitySignal_; }

}
#endif