ddMd/potentials/external/ExternalPotentialImpl.h

#ifndef DDMD_EXTERNAL_POTENTIAL_IMPL_H
#define DDMD_EXTERNAL_POTENTIAL_IMPL_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 <ddMd/potentials/external/ExternalPotential.h>
#include <util/global.h>

namespace Util {
   class Vector;
   class Tensor;
}

namespace DdMd
{

   using namespace Util;

   class Simulation;
   class AtomStorage;

   template <class Interaction>
   class ExternalPotentialImpl : public ExternalPotential
   {

   public:

      ExternalPotentialImpl(Simulation& simulation);

      ExternalPotentialImpl();

      virtual ~ExternalPotentialImpl();

      virtual void associate(Boundary& boundary, AtomStorage& storage);

      virtual void setNAtomType(int nAtomType);

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

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

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



      virtual double externalEnergy(const Vector& position, int typeId) const;

      virtual void getExternalForce(const Vector& position, int typeId, 
                                    Vector& force) const; 

      void set(std::string name, double value)
      {  interactionPtr_->set(name, value); }

      double get(std::string name) const
      {  return interactionPtr_->get(name); }

      virtual std::string interactionClassName() const;

      Interaction& interaction();

      const Interaction& interaction() const;



      virtual void computeForces();

      #ifdef UTIL_MPI
      virtual void computeEnergy(MPI::Intracomm& communicator);
      #else
      virtual void computeEnergy();
      #endif


   private:

      Interaction* interactionPtr_;

      bool isInitialized_;

      double computeForces(bool needForce, bool needEnergy);

   };

}

#include <ddMd/simulation/Simulation.h>
#include <ddMd/storage/AtomStorage.h>
#include <ddMd/storage/AtomIterator.h>

#include <util/space/Dimension.h>
#include <util/space/Vector.h>
#include <util/global.h>

#include <fstream>

namespace DdMd
{

   using namespace Util;

   /* 
   * Constructor.
   */
   template <class Interaction>
   ExternalPotentialImpl<Interaction>::ExternalPotentialImpl(Simulation& simulation)
    : ExternalPotential(simulation),
      interactionPtr_(0),
      isInitialized_(false)
   {  
      interactionPtr_ = new Interaction;
      interaction().setBoundary(simulation.boundary());
      setNAtomType(simulation.nAtomType());
   }
 
   /* 
   * Default constructor.
   */
   template <class Interaction>
   ExternalPotentialImpl<Interaction>::ExternalPotentialImpl()
    : ExternalPotential(),
      interactionPtr_(0),
      isInitialized_(false)
   {  interactionPtr_ = new Interaction; }
 
   /* 
   * Destructor. 
   */
   template <class Interaction>
   ExternalPotentialImpl<Interaction>::~ExternalPotentialImpl() 
   {
      if (interactionPtr_) {
         delete interactionPtr_;
      }
   }

   /*
   * Associate with related objects. (for unit testing).
   *
   * Required iff instantiated with default constructor.
   */
   template <class Interaction>
   void ExternalPotentialImpl<Interaction>
        ::associate(Boundary& boundary, AtomStorage& storage)
   {
      ExternalPotential::associate(boundary, storage);
      interaction().setBoundary(boundary);
   } 

   template <class Interaction>
   void ExternalPotentialImpl<Interaction>::setNAtomType(int nAtomType)
   {  interaction().setNAtomType(nAtomType); }

   /* 
   * Read parameters from file.
   */
   template <class Interaction>
   void ExternalPotentialImpl<Interaction>::readParameters(std::istream &in) 
   {
      UTIL_CHECK(!isInitialized_);
      bool nextIndent = false; // Do not indent interaction block
      addParamComposite(interaction(), nextIndent);
      interaction().readParameters(in);
      isInitialized_ = true;
   }

   /*
   * Load internal state from an archive.
   */
   template <class Interaction>
   void ExternalPotentialImpl<Interaction>::loadParameters(Serializable::IArchive &ar)
   {
      UTIL_CHECK(!isInitialized_);
      bool nextIndent = false; // Do not indent interaction block
      addParamComposite(interaction(), nextIndent);
      interaction().loadParameters(ar);
      isInitialized_ = true;
   }

   /*
   * Save internal state to an archive (call only on master processor).
   */
   template <class Interaction>
   void ExternalPotentialImpl<Interaction>::save(Serializable::OArchive &ar)
   {  interaction().save(ar); }
 
   /*
   * Return external energy of an atom.
   */
   template <class Interaction>
   double 
   ExternalPotentialImpl<Interaction>::externalEnergy(const Vector& position, 
                                                      int typeId) const
   {  return interaction().energy(position, typeId); }

   /*
   * Return external force on an atom.
   */
   template <class Interaction>
   void 
    ExternalPotentialImpl<Interaction>::getExternalForce(const Vector& position, 
                                                int typeId, Vector& force) const
   { interaction().getForce(position, typeId, force); }

   /*
   * Return external interaction class name.
   */
   template <class Interaction>
   std::string ExternalPotentialImpl<Interaction>::interactionClassName() const
   {  return interaction().className(); }

   template <class Interaction>
   const Interaction& ExternalPotentialImpl<Interaction>::interaction() const
   {
      assert(interactionPtr_);  
      return *interactionPtr_; 
   }

   template <class Interaction>
   Interaction& ExternalPotentialImpl<Interaction>::interaction()
   {  
      assert(interactionPtr_);  
      return *interactionPtr_; 
   }

   /*
   * Increment atomic forces, without calculating energy.
   */
   template <class Interaction>
   void ExternalPotentialImpl<Interaction>::computeForces()
   {  computeForces(true, false); }

   /*
   * Compute total external energy on all processors.
   */
   template <class Interaction>
   #ifdef UTIL_MPI
   void 
   ExternalPotentialImpl<Interaction>::computeEnergy(MPI::Intracomm& communicator)
   #else
   void ExternalPotentialImpl<Interaction>::computeEnergy()
   #endif
   { 
      double localEnergy = 0; 
      localEnergy = computeForces(false, true); 
      #ifdef UTIL_MPI
      reduceEnergy(localEnergy, communicator);
      #else
      setEnergy(localEnergy);
      #endif
   }

   /*
   * Increment atomic forces and/or external energy (private).
   */
   template <class Interaction>
   double 
   ExternalPotentialImpl<Interaction>::computeForces(bool needForce, bool needEnergy)
   {
      Vector f;
      double energy = 0.0;
      AtomIterator iter;
      int type;

      storage().begin(iter);
      for ( ; iter.notEnd(); ++iter) {
         type = iter->typeId();
         if (needEnergy) {
            energy += 
               interactionPtr_->energy(iter->position(), type);
         }
         if (needForce) {
            interactionPtr_->getForce(iter->position(), type, f);
            iter->force() += f;
         }
      }
      return energy;
   }

}
#endif