ddMd/potentials/angle/AnglePotentialImpl.h

#ifndef DDMD_ANGLE_POTENTIAL_IMPL_H
#define DDMD_ANGLE_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 "AnglePotential.h" // base class

namespace Util
{
   class Vector;
   class Tensor;
}

namespace DdMd
{

   using namespace Util;

   class Simulation;
   template <int N> class GroupStorage;

   template <class Interaction>
   class AnglePotentialImpl : public AnglePotential
   {

   public:

      AnglePotentialImpl(Simulation& simulation);

      AnglePotentialImpl();

      virtual ~AnglePotentialImpl();

      virtual void setNAngleType(int nAtomType);
  
      virtual void readParameters(std::istream& in);

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

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



      double angleEnergy(double cosTheta, int type) const;
 
      void angleForce(const Vector& R1, const Vector& R2,
                      Vector& F1, Vector& F2, int type) const;

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

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

      virtual std::string interactionClassName() const;

      const Interaction& interaction() const;

      Interaction& interaction();



      virtual void computeForces();

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

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


   private:

      Interaction* interactionPtr_;

      bool isInitialized_;

   };

}

#include "AnglePotential.h"
#include <ddMd/simulation/Simulation.h>
#include <ddMd/storage/GroupStorage.h>
#include <ddMd/storage/GroupIterator.h>

#include <simp/boundary/Boundary.h>
#include <util/space/Dimension.h>
#include <util/space/Vector.h>
#include <util/space/Tensor.h>
#include <util/global.h>

#include <fstream>

namespace DdMd
{

   using namespace Util;
   using namespace Simp;

   /* 
   * Constructor.
   */
   template <class Interaction>
   AnglePotentialImpl<Interaction>::AnglePotentialImpl(Simulation& simulation)
    : AnglePotential(simulation),
      interactionPtr_(0),
      isInitialized_(false)
   {  
      interactionPtr_ = new Interaction(); 
      setNAngleType(simulation.nAngleType());
   }
 
   /* 
   * Default constructor.
   */
   template <class Interaction>
   AnglePotentialImpl<Interaction>::AnglePotentialImpl()
    : AnglePotential(),
      interactionPtr_(0),
      isInitialized_(false)
   {  interactionPtr_ = new Interaction(); }
 
   /* 
   * Destructor. 
   */
   template <class Interaction>
   AnglePotentialImpl<Interaction>::~AnglePotentialImpl() 
   {
      if (interactionPtr_) {
         delete interactionPtr_;
      }
   }

   /*
   * Set the maximum number of atom types.
   */
   template <class Interaction>
   void AnglePotentialImpl<Interaction>::setNAngleType(int nAngleType)
   {  interaction().setNAngleType(nAngleType); }

   /*
   * Read angle interaction parameters.
   */
   template <class Interaction>
   void AnglePotentialImpl<Interaction>::readParameters(std::istream& in)
   {
      UTIL_CHECK(!isInitialized_);
      bool nextIndent = false;
      addParamComposite(interaction(), nextIndent);
      interaction().readParameters(in);
      isInitialized_ = true;
   }

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

   /*
   * Save internal state to an archive.
   */
   template <class Interaction>
   void AnglePotentialImpl<Interaction>::save(Serializable::OArchive &ar)
   {  interaction().save(ar); }

   /*
   * Return energy for a single angle.
   */
   template <class Interaction>
   double 
   AnglePotentialImpl<Interaction>::angleEnergy(double cosTheta, int angleTypeId) 
      const
   {  return interaction().energy(cosTheta, angleTypeId); }

   /*
   * Return forces for a single angle.
   */
   template <class Interaction> void 
   AnglePotentialImpl<Interaction>::angleForce(const Vector& R1, 
                                               const Vector& R2, 
                                               Vector& F1, Vector& F2, 
                                               int typeId) const
   {  interaction().force(R1, R2, F1, F2, typeId); }

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

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

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

   /*
   * Increment atomic forces, without calculating energy.
   */
   template <class Interaction>
   void AnglePotentialImpl<Interaction>::computeForces()
   {  
      Vector dr1; // R[1] - R[0]
      Vector dr2; // R[2] - R[1]
      Vector f1;  // d(energy)/d(dr1)
      Vector f2;  // d(energy)/d(dr2)
      GroupIterator<3> iter;
      Atom* atom0Ptr;
      Atom* atom1Ptr;
      Atom* atom2Ptr;
      int   type;

      storage().begin(iter);
      for ( ; iter.notEnd(); ++iter) {
         type = iter->typeId();
         atom0Ptr = iter->atomPtr(0);
         atom1Ptr = iter->atomPtr(1);
         atom2Ptr = iter->atomPtr(2);
         // Calculate minimimum image separations
         boundary().distanceSq(atom1Ptr->position(),
                               atom0Ptr->position(), dr1);
         boundary().distanceSq(atom2Ptr->position(),
                               atom1Ptr->position(), dr2);
         interaction().force(dr1, dr2, f1, f2, type);
         if (!atom0Ptr->isGhost()) {
            atom0Ptr->force() += f1;
         }
         if (!atom1Ptr->isGhost()) {
            atom1Ptr->force() -= f1;
            atom1Ptr->force() += f2;
         }
         if (!atom2Ptr->isGhost()) {
            atom2Ptr->force() -= f2;
         }
      }
   }

   /*
   * Compute total angle energy on all processors.
   */
   template <class Interaction>
   #ifdef UTIL_MPI
   void 
   AnglePotentialImpl<Interaction>::computeEnergy(MPI::Intracomm& communicator)
   #else
   void AnglePotentialImpl<Interaction>::computeEnergy()
   #endif
   {

      // Do nothing and return if energy is already set.
      if (isEnergySet()) return;
 
      Vector dr1; // R[1] - R[0]
      Vector dr2; // R[2] - R[1]
      double rsq1, rsq2, cosTheta;
      double angleEnergy;
      double localEnergy = 0.0;
      double fraction;
      double third = 1.0/3.0;
      GroupIterator<3> iter;
      Atom* atom0Ptr;
      Atom* atom1Ptr;
      Atom* atom2Ptr;
      int   type, isLocal0, isLocal1, isLocal2;

      storage().begin(iter);
      for ( ; iter.notEnd(); ++iter) {
         type = iter->typeId();
         atom0Ptr = iter->atomPtr(0);
         atom1Ptr = iter->atomPtr(1);
         atom2Ptr = iter->atomPtr(2);
         isLocal0 = !(atom0Ptr->isGhost());
         isLocal1 = !(atom1Ptr->isGhost());
         isLocal2 = !(atom2Ptr->isGhost());
         // Calculate minimimum image separations
         rsq1 = boundary().distanceSq(atom1Ptr->position(),
                                      atom0Ptr->position(), dr1);
         rsq2 = boundary().distanceSq(atom2Ptr->position(),
                                      atom1Ptr->position(), dr2);
         cosTheta = dr1.dot(dr2) / sqrt(rsq1 * rsq2);
         angleEnergy = interaction().energy(cosTheta, type);
         fraction = (isLocal0 + isLocal1 + isLocal2)*third;
         localEnergy += fraction*angleEnergy;
      }

      // Add localEnergy from all nodes, set sum on master.
      reduceEnergy(localEnergy, communicator);
   }

   /*
   * Compute total pair stress (Call on all processors).
   */
   template <class Interaction>
   #ifdef UTIL_MPI
   void AnglePotentialImpl<Interaction>::computeStress(MPI::Intracomm& communicator)
   #else
   void AnglePotentialImpl<Interaction>::computeStress()
   #endif
   {
      // Do nothing and return if stress is already set.
      if (isStressSet()) return;
 
      Tensor localStress;
      Vector dr1, dr2;
      Vector f1, f2;
      double factor;
      double prefactor = -1.0/3.0;
      GroupIterator<3> iter;
      Atom*  atom0Ptr;
      Atom*  atom1Ptr;
      Atom*  atom2Ptr;
      int    type;
      int    isLocal0, isLocal1, isLocal2;

      // Iterate over angle groups
      localStress.zero();
      storage().begin(iter);
      for ( ; iter.notEnd(); ++iter) {
         atom0Ptr = iter->atomPtr(0);
         atom1Ptr = iter->atomPtr(1);
         atom2Ptr = iter->atomPtr(2);
         type = iter->typeId();
         boundary().distanceSq(atom1Ptr->position(),
                                  atom0Ptr->position(), dr1);
         boundary().distanceSq(atom2Ptr->position(),
                                  atom1Ptr->position(), dr2);

         // Calculate derivatives f1, f2 of energy with respect to dr1, dr2
         interaction().force(dr1, dr2, f1, f2, type);

         isLocal0 = !(atom0Ptr->isGhost());
         isLocal1 = !(atom1Ptr->isGhost());
         isLocal2 = !(atom2Ptr->isGhost());
         factor = prefactor*(isLocal0 + isLocal1 + isLocal2);

         // Increment localStress tensor
         dr1 *= factor;
         dr2 *= factor;
         incrementPairStress(f1, dr1, localStress);
         incrementPairStress(f2, dr2, localStress);
      }

      // Normalize by volume 
      localStress /= boundary().volume();

      // Add localEnergy from all nodes, set sum on master.
      reduceStress(localStress, communicator);
   }

}
#endif