Line data    Source code

       1              : // Copyright (C) 2015 Technische Universitaet Muenchen
       2              : // This file is part of the Mamico project. For conditions of distribution
       3              : // and use, please see the copyright notice in Mamico's main folder, or at
       4              : // www5.in.tum.de/mamico
       5              : #ifndef _MOLECULARDYNAMICS_COUPLING_KINETICENERGYCONTROLLER_H_
       6              : #define _MOLECULARDYNAMICS_COUPLING_KINETICENERGYCONTROLLER_H_
       7              : 
       8              : #include "coupling/cell-mappings/ComputeKineticEnergyMapping.h"
       9              : #include "coupling/cell-mappings/ComputeMassMapping.h"
      10              : #include "coupling/cell-mappings/ComputeMomentumMapping.h"
      11              : #include "coupling/cell-mappings/ComputeTemperatureMapping.h"
      12              : #include "coupling/cell-mappings/SetKineticEnergyMapping.h"
      13              : #include "coupling/cell-mappings/SetTemperatureMapping.h"
      14              : #include "coupling/datastructures/CouplingCell.h"
      15              : #include "coupling/datastructures/CouplingCellWithLinkedCells.h"
      16              : 
      17              : namespace coupling {
      18              : template <class LinkedCell, unsigned int dim> class KineticEnergyController;
      19              : }
      20              : 
      21              : /** @brief controles and regulates the kinetic energy of the MD system.
      22              :  *  @author Philipp Neumann
      23              :  *  @tparam LinkedCell the LinkedCell class is given by the implementation of
      24              :  * linked cells in the molecular dynamics simulation
      25              :  *  @tparam dim  refers to the spacial dimension of the simulation, can be 1, 2,
      26              :  * or 3*/
      27              : template <class LinkedCell, unsigned int dim> class coupling::KineticEnergyController {
      28              : public:
      29              :   /** @brief a simple constructor
      30              :    *  @param mdSolverInterface interface to the md solver */
      31            4 :   KineticEnergyController(coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface) : _mdSolverInterface(mdSolverInterface) {}
      32              :   /** @brief a simple destructor*/
      33            0 :   ~KineticEnergyController() {}
      34              : 
      35              :   /** @brief computes and returns the kinetic energy within a coupling cell
      36              :    *  @param cell coupling cell to compute the kinetic energy for
      37              :    *  @returns the kinetic energy in the cell */
      38              :   double computeKineticEnergy(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell) const {
      39              :     coupling::cellmappings::ComputeKineticEnergyMapping<LinkedCell, dim> computeKineticEnergyMapping(_mdSolverInterface);
      40              :     cell.iterateConstCells(computeKineticEnergyMapping);
      41              :     return computeKineticEnergyMapping.getKineticEnergy();
      42              :   }
      43              : 
      44              :   /** @brief computes and returns the temperature within a coupling cell
      45              :    *  @param cell coupling cell to compute the temperature in
      46              :    *  @returns the temperature in the cell */
      47            0 :   double computeTemperature(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell) const {
      48            0 :     coupling::cellmappings::ComputeMomentumMapping<LinkedCell, dim> computeMomentumMapping(_mdSolverInterface);
      49            0 :     cell.iterateConstCells(computeMomentumMapping);
      50            0 :     tarch::la::Vector<dim, double> meanVelocity = computeMomentumMapping.getMeanVelocity();
      51              : 
      52            0 :     coupling::cellmappings::ComputeTemperatureMapping<LinkedCell, dim> computeTemperatureMapping(meanVelocity, _mdSolverInterface);
      53            0 :     cell.iterateConstCells(computeTemperatureMapping);
      54            0 :     return computeTemperatureMapping.getTemperature();
      55            0 :   }
      56              : 
      57              :   /** Therefore, the mean velocity is computed first. Afterwards the deviation
      58              :    * from the mean velocity are rescaled such that momentum is conserved.
      59              :    *  @brief sets the kinetic energy within a coupling cell to the input
      60              :    * value.
      61              :    *  @param cell the coupling cell to set the kinetic energy in
      62              :    *  @param kineticEnergy the value the kinetic energy shall be set to*/
      63              :   void setKineticEnergy(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const double& kineticEnergy) const {
      64              :     // determine mass, momentum and old kinetic energy
      65              :     coupling::cellmappings::ComputeMassMapping<LinkedCell, dim> computeMassMapping(_mdSolverInterface);
      66              :     cell.iterateConstCells(computeMassMapping);
      67              :     coupling::cellmappings::ComputeMomentumMapping<LinkedCell, dim> computeMomentumMapping(_mdSolverInterface);
      68              :     cell.iterateConstCells(computeMomentumMapping);
      69              :     coupling::cellmappings::ComputeKineticEnergyMapping<LinkedCell, dim> computeKineticEnergyMapping(_mdSolverInterface);
      70              :     cell.iterateConstCells(computeKineticEnergyMapping);
      71              :     // set new kinetic energy
      72              :     unsigned int numberParticles = computeMassMapping.getNumberOfParticles();
      73              :     tarch::la::Vector<dim, double> meanVelocity = computeMomentumMapping.getMeanVelocity();
      74              :     double oldKineticEnergy = computeKineticEnergyMapping.getKineticEnergy();
      75              :     coupling::cellmappings::SetKineticEnergyMapping<LinkedCell, dim> setKineticEnergyMapping(oldKineticEnergy, kineticEnergy, numberParticles, meanVelocity);
      76              :     cell.iterateCells(setKineticEnergyMapping);
      77              :   }
      78              : 
      79              :   /** Here we just scale the deviation from the mean velocity accordingly, i.e.
      80              :    * we set: v_molecule = v_mean +
      81              :    * sqrt(temperature/current_temperature)*(v_molecule-v_mean)
      82              :    *  @brief sets the temperature within the cell to 'temperature'.
      83              :    *  @param cell the coupling cell the temperature shall be applied in
      84              :    *  @param temperature the value the temperature shall be set at*/
      85            0 :   void setTemperature(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const double& temperature) const {
      86            0 :     coupling::cellmappings::ComputeMomentumMapping<LinkedCell, dim> computeMomentumMapping(_mdSolverInterface);
      87            0 :     cell.iterateConstCells(computeMomentumMapping);
      88            0 :     tarch::la::Vector<dim, double> meanVelocity = computeMomentumMapping.getMeanVelocity();
      89              : 
      90            0 :     coupling::cellmappings::ComputeTemperatureMapping<LinkedCell, dim> computeTemperatureMapping(meanVelocity, _mdSolverInterface);
      91            0 :     cell.iterateConstCells(computeTemperatureMapping);
      92            0 :     double currentTemperature = computeTemperatureMapping.getTemperature();
      93              : 
      94            0 :     coupling::cellmappings::SetTemperatureMapping<LinkedCell, dim> setTemperatureMapping(currentTemperature, temperature, meanVelocity, _mdSolverInterface);
      95            0 :     cell.iterateCells(setTemperatureMapping);
      96            0 :   }
      97              : 
      98              : private:
      99              :   /** interface of the md solver*/
     100              :   coupling::interface::MDSolverInterface<LinkedCell, dim>* const _mdSolverInterface;
     101              : };
     102              : #endif // _MOLECULARDYNAMICS_COUPLING_KINETICENERGYCONTROLLER_H_