ddMd/potentials/dihedral/DihedralPotentialImpl.h

#ifndef DDMD_DIHEDRAL_POTENTIAL_IMPL_H
#define DDMD_DIHEDRAL_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 <util/global.h>
#include "DihedralPotential.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 DihedralPotentialImpl : public DihedralPotential
   {

   public:

      DihedralPotentialImpl(Simulation& simulation);

      DihedralPotentialImpl();

      virtual ~DihedralPotentialImpl();



      virtual void setNDihedralType(int nAtomType);

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

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

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



      void set(std::string name, int type, double value);

      double get(std::string name, int type) const;

      virtual double
      dihedralEnergy(const Vector& R1, const Vector& R2, const Vector& R3,
                     int type) const;

      virtual void
      dihedralForce(const Vector& R1, const Vector& R2, const Vector& R3,
                    Vector& F1, Vector& F2, Vector& F3, int type) const;

      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

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


   private:

      Interaction* interactionPtr_;

      bool isInitialized_;

   };

}

#include "DihedralPotential.h"
#include <ddMd/simulation/Simulation.h>
#include <ddMd/storage/GroupStorage.h>
#include <ddMd/storage/GroupIterator.h>
#include <simp/boundary/Boundary.h>
#include <util/space/Vector.h>
#include <util/space/Dimension.h>
#include <util/space/Vector.h>
#include <util/space/Tensor.h>
//#include <util/accumulators/setToZero.h>
#include <util/global.h>

#include <fstream>

namespace DdMd
{

   using namespace Util;
   using namespace Simp;

   /*
   * Constructor.
   */
   template <class Interaction>
   DihedralPotentialImpl<Interaction>::DihedralPotentialImpl(Simulation& simulation)
    : DihedralPotential(simulation),
      interactionPtr_(0),
      isInitialized_(false)
   {
      interactionPtr_ = new Interaction();
      setNDihedralType(simulation.nDihedralType());
   }

   /*
   * Default constructor.
   */
   template <class Interaction>
   DihedralPotentialImpl<Interaction>::DihedralPotentialImpl()
    : DihedralPotential(),
      interactionPtr_(0),
      isInitialized_(false)
   {  interactionPtr_ = new Interaction(); }

   /*
   * Destructor.
   */
   template <class Interaction>
   DihedralPotentialImpl<Interaction>::~DihedralPotentialImpl()
   {
      if (interactionPtr_) {
         delete interactionPtr_;
      }
   }

   /*
   * Set the maximum number of dihedral types.
   */
   template <class Interaction>
   void DihedralPotentialImpl<Interaction>::setNDihedralType(int nDihedralType)
   {  interaction().setNDihedralType(nDihedralType); }

   /*
   * Read dihedral interaction parameters.
   */
   template <class Interaction>
   void DihedralPotentialImpl<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
   DihedralPotentialImpl<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 DihedralPotentialImpl<Interaction>::save(Serializable::OArchive &ar)
   {  interaction().save(ar); }

   /*
   * Modify an dihedral parameter, identified by a string.
   */
   template <class Interaction>
   inline void DihedralPotentialImpl<Interaction>::
      set(std::string name, int type, double value)
   {  interactionPtr_->set(name, type, value); }

   /*
   * Get a parameter value, identified by a string.
   */
   template <class Interaction>
   inline double DihedralPotentialImpl<Interaction>::
      get(std::string name, int type) const
   {  return interactionPtr_->get(name, type); }

   /*
   * Return energy for a single dihedral.
   */
   template <class Interaction>
   inline double DihedralPotentialImpl<Interaction>::
      dihedralEnergy(const Vector& R1, const Vector& R2, const Vector& R3,
                     int typeId) const
   {  return interaction().energy(R1, R2, R3, typeId); }

   /*
   * Return forces for a single dihedral.
   */
   template <class Interaction>
   inline void DihedralPotentialImpl<Interaction>::
      dihedralForce(const Vector& R1, const Vector& R2, const Vector& R3,
                    Vector& F1, Vector& F2, Vector& F3, int typeId) const
   {  interaction().force(R1, R2, R3, F1, F2, F3, typeId); }

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

   /*
   * Get Interaction by reference.
   */
   template <class Interaction>
   inline Interaction& DihedralPotentialImpl<Interaction>::interaction()
   {  return *interactionPtr_; }

   /*
   * Get Interaction by const reference.
   */
   template <class Interaction>
   inline const Interaction& DihedralPotentialImpl<Interaction>::interaction() const
   {  return *interactionPtr_; }

   /*
   * Increment atomic forces, without calculating energy.
   */
   template <class Interaction>
   void DihedralPotentialImpl<Interaction>::computeForces()
   {
      UTIL_CHECK(isInitialized_);

      Vector dr1; // R[1] - R[0]
      Vector dr2; // R[2] - R[1]
      Vector dr3; // R[3] - R[2]
      Vector f1, f2, f3;
      GroupIterator<4> iter;
      Atom* atom0Ptr;
      Atom* atom1Ptr;
      Atom* atom2Ptr;
      Atom* atom3Ptr;
      int   type;

      // Loop over dihedral groups
      for (storage().begin(iter); iter.notEnd(); ++iter) {
         type = iter->typeId();
         atom0Ptr = iter->atomPtr(0);
         atom1Ptr = iter->atomPtr(1);
         atom2Ptr = iter->atomPtr(2);
         atom3Ptr = iter->atomPtr(3);
         // Calculate minimimum image separations dr1, dr2, dr3
         boundary().distanceSq(atom1Ptr->position(),
                               atom0Ptr->position(), dr1);
         boundary().distanceSq(atom2Ptr->position(),
                               atom1Ptr->position(), dr2);
         boundary().distanceSq(atom3Ptr->position(),
                               atom2Ptr->position(), dr3);

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

         if (!atom0Ptr->isGhost()) {
            atom0Ptr->force() += f1;
         }
         if (!atom1Ptr->isGhost()) {
            atom1Ptr->force() -= f1;
            atom1Ptr->force() += f2;
         }
         if (!atom2Ptr->isGhost()) {
            atom2Ptr->force() -= f2;
            atom2Ptr->force() += f3;
         }
         if (!atom3Ptr->isGhost()) {
            atom3Ptr->force() -= f3;
         }
      }
   }

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

      // If energy is already set, do nothing and return.
      if (isEnergySet()) return;

      Vector dr1; // R[1] - R[0]
      Vector dr2; // R[2] - R[1]
      Vector dr3; // R[3] - R[2]
      double dihedralEnergy;
      double localEnergy;
      double fraction;
      GroupIterator<4> iter;
      Atom* atom0Ptr;
      Atom* atom1Ptr;
      Atom* atom2Ptr;
      Atom* atom3Ptr;
      int   type, isLocal0, isLocal1, isLocal2, isLocal3;

      // Loop over dihedral groups
      localEnergy = 0.0;
      for (storage().begin(iter) ; iter.notEnd(); ++iter) {
         type = iter->typeId();
         atom0Ptr = iter->atomPtr(0);
         atom1Ptr = iter->atomPtr(1);
         atom2Ptr = iter->atomPtr(2);
         atom3Ptr = iter->atomPtr(3);

         // Calculate minimimum image separation vectors
         boundary().distanceSq(atom1Ptr->position(),
                               atom0Ptr->position(), dr1);
         boundary().distanceSq(atom2Ptr->position(),
                               atom1Ptr->position(), dr2);
         boundary().distanceSq(atom3Ptr->position(),
                               atom2Ptr->position(), dr3);

         // Calculate energy for one dihedral
         dihedralEnergy = interaction().energy(dr1, dr2, dr3, type);

         isLocal0 = !(atom0Ptr->isGhost());
         isLocal1 = !(atom1Ptr->isGhost());
         isLocal2 = !(atom2Ptr->isGhost());
         isLocal3 = !(atom3Ptr->isGhost());
         fraction = (isLocal0 + isLocal1 + isLocal2 + isLocal3)*0.25;
         localEnergy += fraction*dihedralEnergy;
      }

      // 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
   DihedralPotentialImpl<Interaction>::computeStress(MPI::Intracomm& communicator)
   #else
   void DihedralPotentialImpl<Interaction>::computeStress()
   #endif
   {
      // If stress is already set, do nothing and return.
      if (isStressSet()) return;

      Tensor localStress;
      Vector dr1, dr2, dr3;
      Vector f1,  f2, f3;
      double factor;
      double prefactor = -0.25;
      GroupIterator<4> iter;
      Atom*  atom0Ptr;
      Atom*  atom1Ptr;
      Atom*  atom2Ptr;
      Atom*  atom3Ptr;
      int    type;
      int    isLocal0, isLocal1, isLocal2, isLocal3;

      // Iterate over dihedrals
      localStress.zero();
      for (storage().begin(iter); iter.notEnd(); ++iter) {
         type = iter->typeId();
         atom0Ptr = iter->atomPtr(0);
         atom1Ptr = iter->atomPtr(1);
         atom2Ptr = iter->atomPtr(2);
         atom3Ptr = iter->atomPtr(3);

         // Calculate stress here
         boundary().distanceSq(atom1Ptr->position(),
                               atom0Ptr->position(), dr1);
         boundary().distanceSq(atom2Ptr->position(),
                               atom1Ptr->position(), dr2);
         boundary().distanceSq(atom3Ptr->position(),
                               atom2Ptr->position(), dr3);

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

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

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

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

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

}
#endif