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_INTERFACE_SIMPLEMDSOLVERINTERFACE_CPP_
       6             : #define _MOLECULARDYNAMICS_COUPLING_INTERFACE_SIMPLEMDSOLVERINTERFACE_CPP_
       7             : #include "simplemd/BoundaryTreatment.h"
       8             : #include "simplemd/MolecularDynamicsDefinitions.h"
       9             : #include "simplemd/cell-mappings/LennardJonesForceMapping.h"
      10             : #include "simplemd/cell-mappings/LennardJonesPotentialEnergyMapping.h"
      11             : #include "simplemd/cell-mappings/ResetPotentialEnergyMapping.h"
      12             : #include "simplemd/services/ExternalForceService.h"
      13             : #include "simplemd/services/LinkedCellService.h"
      14             : #include "simplemd/services/MolecularPropertiesService.h"
      15             : #include "simplemd/services/ParallelTopologyService.h"
      16             : 
      17             : #include "coupling/CouplingMDDefinitions.h"
      18             : #include "coupling/interface/MDSolverInterface.h"
      19             : #include "coupling/interface/impl/SimpleMD/SimpleMDMoleculeIterator.h"
      20             : 
      21             : namespace coupling {
      22             : namespace interface {
      23             : 
      24             : /** general MD solver interface for SimpleMD.
      25             :  *  @author Philipp Neumann
      26             :  */
      27             : class SimpleMDSolverInterface : public MDSolverInterface<simplemd::LinkedCell, MD_DIM> {
      28             : private:
      29             :   simplemd::services::ParallelTopologyService& _parallelTopologyService;
      30             :   simplemd::services::MoleculeService& _moleculeService;
      31             :   simplemd::services::LinkedCellService& _linkedCellService;
      32             :   const simplemd::services::MolecularPropertiesService& _molecularPropertiesService;
      33             :   /** used for the synchronization of molecules in boundary regions and between different processes when
      34             :    *  mass was inserted/ deleted by the coupling.
      35             :    */
      36             :   simplemd::BoundaryTreatment& _boundaryTreatment;
      37             :   const tarch::la::Vector<MD_LINKED_CELL_NEIGHBOURS, simplemd::BoundaryType>& _localBoundaryInformation;
      38             : 
      39             :   /** number of linked cells in each coupling cell */
      40             : 
      41             :   const double _sigma6;
      42             :   const double _epsilon;
      43             :   const double _cutoffRadiusSquared;
      44             :   const double _cutoffEnergy;
      45             : 
      46             :   const double _dt;
      47             : 
      48             :   tarch::la::Vector<MD_DIM, unsigned int>
      49             :   computeNumberOfLinkedCellsPerCouplingCells(const tarch::la::Vector<MD_DIM, unsigned int>& localNumberOfCouplingCells) const {
      50             :     tarch::la::Vector<MD_DIM, unsigned int> numberOfLinkedCells(0);
      51             :     for (unsigned int d = 0; d < MD_DIM; d++) {
      52             :       numberOfLinkedCells[d] = _linkedCellService.getLocalNumberOfCells()[d] / localNumberOfCouplingCells[d];
      53             : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
      54             :       if (numberOfLinkedCells[d] * localNumberOfCouplingCells[d] != _linkedCellService.getLocalNumberOfCells()[d]) {
      55             :         std::cout << "ERROR coupling::interface::SimpleMDSolverInterface: Number of linked cells is not a multiple of the coupling cells!" << std::endl;
      56             :         std::cout << "Number linked cells: " << numberOfLinkedCells << ", coupling cells: " << localNumberOfCouplingCells << std::endl;
      57             :         std::cout << "Local number of linked cells: " << _linkedCellService.getLocalNumberOfCells() << std::endl;
      58             :         exit(EXIT_FAILURE);
      59             :       }
      60             : #endif
      61             :     }
      62             :     return numberOfLinkedCells;
      63             :   }
      64             : 
      65             :   /** Though this method is available in LennardJonesPotentialEnergyMapping, I re-implemented it at
      66             :    *  this stage, as the latter always fetches the epsilon,sigma-parameters first. This is nice for
      67             :    *  simulations where parameters are changed over time; however for static simulations it just slows
      68             :    *  down...
      69             :    */
      70           0 :   double getPotentialEnergy(const tarch::la::Vector<MD_DIM, double>& position1, const tarch::la::Vector<MD_DIM, double>& position2) const {
      71           0 :     const tarch::la::Vector<MD_DIM, double> distance(position2 - position1);
      72           0 :     const double r2 = tarch::la::dot(distance, distance);
      73             :     // if we are in the near distance, do computation, otherwise return 0.0
      74           0 :     if (r2 <= _cutoffRadiusSquared) {
      75           0 :       const double r6 = r2 * r2 * r2;
      76           0 :       const double energyBuffer = 4.0 * _epsilon * (_sigma6 / r6) * ((_sigma6 / r6) - 1.0) - _cutoffEnergy;
      77           0 :       return 0.5 * energyBuffer;
      78             :     } else {
      79             :       return 0.0;
      80             :     }
      81             :   }
      82             : 
      83           0 :   tarch::la::Vector<MD_DIM, double> getLennardJonesForce(const tarch::la::Vector<MD_DIM, double>& position1,
      84             :                                                          const tarch::la::Vector<MD_DIM, double>& position2) const {
      85           0 :     const tarch::la::Vector<MD_DIM, double> rij(position2 - position1);
      86           0 :     const double rij2 = tarch::la::dot(rij, rij);
      87             : #if (MD_ERROR == MD_YES)
      88             :     if (rij2 == 0.0) {
      89             :       std::cout << "ERROR cellmappings::LennardJonesForceMapping::getLennardJonesForce(): Particle positions are identical!" << std::endl;
      90             :       exit(EXIT_FAILURE);
      91             :     }
      92             : #endif
      93             : 
      94           0 :     if (rij2 <= _cutoffRadiusSquared) {
      95           0 :       const double rij6 = rij2 * rij2 * rij2;
      96           0 :       return 24.0 * _epsilon / rij2 * (_sigma6 / rij6) * (1.0 - 2.0 * (_sigma6 / rij6)) * rij;
      97             :     } else {
      98           0 :       return tarch::la::Vector<MD_DIM, double>(0.0);
      99             :     }
     100             :   }
     101             : 
     102             : public:
     103          12 :   SimpleMDSolverInterface(simplemd::BoundaryTreatment& boundaryTreatment,
     104             :                           simplemd::services::ParallelTopologyService& parallelTopologyService, simplemd::services::MoleculeService& moleculeService,
     105             :                           simplemd::services::LinkedCellService& linkedCellService,
     106             :                           const simplemd::services::MolecularPropertiesService& molecularPropertiesService,
     107             :                           const tarch::la::Vector<MD_LINKED_CELL_NEIGHBOURS, simplemd::BoundaryType>& localBoundaryInformation, const double& dt)
     108          12 :       : _parallelTopologyService(parallelTopologyService), _moleculeService(moleculeService), _linkedCellService(linkedCellService),
     109          12 :         _molecularPropertiesService(molecularPropertiesService), _boundaryTreatment(boundaryTreatment), _localBoundaryInformation(localBoundaryInformation),
     110          12 :         _sigma6(molecularPropertiesService.getMolecularProperties().getSigma() * molecularPropertiesService.getMolecularProperties().getSigma() *
     111          12 :                 molecularPropertiesService.getMolecularProperties().getSigma() * molecularPropertiesService.getMolecularProperties().getSigma() *
     112          12 :                 molecularPropertiesService.getMolecularProperties().getSigma() * molecularPropertiesService.getMolecularProperties().getSigma()),
     113          12 :         _epsilon(molecularPropertiesService.getMolecularProperties().getEpsilon()),
     114          12 :         _cutoffRadiusSquared(molecularPropertiesService.getMolecularProperties().getCutOffRadius() *
     115          12 :                              molecularPropertiesService.getMolecularProperties().getCutOffRadius()),
     116          12 :         _cutoffEnergy(4.0 * _epsilon * _sigma6 / (_cutoffRadiusSquared * _cutoffRadiusSquared * _cutoffRadiusSquared) *
     117          12 :                       (_sigma6 / (_cutoffRadiusSquared * _cutoffRadiusSquared * _cutoffRadiusSquared) - 1.0)),
     118          24 :         _dt(dt) {}
     119          12 :   ~SimpleMDSolverInterface() {}
     120             : 
     121           0 :   simplemd::LinkedCell& getLinkedCell(const CellIndex_T& couplingCellIndex, const tarch::la::Vector<MD_DIM, unsigned int>& linkedCellInCouplingCell,
     122             :                                       const tarch::la::Vector<MD_DIM, unsigned int>& linkedCellsPerCouplingCell) {
     123             :     // no linked cells found in outer region!
     124             : 
     125           0 :     tarch::la::Vector<MD_DIM, unsigned int> index(_linkedCellService.getLocalIndexOfFirstCell());
     126           0 :     for (unsigned int d = 0; d < MD_DIM; d++) {
     127           0 :       index[d] = index[d] + couplingCellIndex.get()[d] * linkedCellsPerCouplingCell[d] + linkedCellInCouplingCell[d];
     128             :     }
     129           0 :     return _linkedCellService.getLinkedCell(index);
     130             :   }
     131             : 
     132             :   /** returns the global size of the box-shaped MD domain */
     133           0 :   tarch::la::Vector<MD_DIM, double> getGlobalMDDomainSize() const { return _parallelTopologyService.getGlobalDomainSize(); }
     134             : 
     135             :   /** returns the offset (i.e. lower,left corner) of MD domain */
     136           0 :   tarch::la::Vector<MD_DIM, double> getGlobalMDDomainOffset() const { return _parallelTopologyService.getGlobalDomainOffset(); }
     137             : 
     138           0 :   double getMoleculeMass() const { return _molecularPropertiesService.getMolecularProperties().getMass(); }
     139             : 
     140           0 :   double getKB() const { return _molecularPropertiesService.getMolecularProperties().getKB(); }
     141             : 
     142           0 :   double getMoleculeSigma() const { return _molecularPropertiesService.getMolecularProperties().getSigma(); }
     143             : 
     144           0 :   double getMoleculeEpsilon() const { return _molecularPropertiesService.getMolecularProperties().getEpsilon(); }
     145             : 
     146           0 :   void getInitialVelocity(const tarch::la::Vector<MD_DIM, double>& meanVelocity, const double& kB, const double& temperature,
     147             :                           tarch::la::Vector<MD_DIM, double>& initialVelocity) const {
     148           0 :     _moleculeService.getInitialVelocity(meanVelocity, kB, temperature, _molecularPropertiesService, initialVelocity);
     149           0 :   }
     150             : 
     151           0 :   void deleteMoleculeFromMDSimulation(const coupling::interface::Molecule<MD_DIM>& molecule, simplemd::LinkedCell& cell) {
     152           0 :     std::list<simplemd::Molecule*>::iterator it = cell.begin();
     153           0 :     const tarch::la::Vector<MD_DIM, double> moleculePosition = molecule.getPosition();
     154           0 :     while (it != cell.end()) {
     155           0 :       const tarch::la::Vector<MD_DIM, double>& itPosition((*it)->getConstPosition());
     156           0 :       if (moleculePosition == itPosition) {
     157           0 :         simplemd::Molecule* myMolecule = (*it);
     158           0 :         cell.deleteMolecule(myMolecule);
     159           0 :         _moleculeService.deleteMolecule(*myMolecule);
     160           0 :         return;
     161             :       }
     162           0 :       it++;
     163             :     }
     164             : 
     165           0 :     std::cout << "Could delete molecule at position " << moleculePosition << "!" << std::endl;
     166           0 :     exit(EXIT_FAILURE);
     167             :   }
     168             : 
     169           0 :   void addMoleculeToMDSimulation(const coupling::interface::Molecule<MD_DIM>& molecule) {
     170           0 :     const tarch::la::Vector<MD_DIM, double> position = molecule.getPosition();
     171           0 :     const tarch::la::Vector<MD_DIM, double> velocity = molecule.getVelocity();
     172           0 :     const tarch::la::Vector<MD_DIM, double> force = molecule.getForce();
     173           0 :     const double potentialEnergy = molecule.getPotentialEnergy();
     174           0 :     simplemd::Molecule newMolecule(position, velocity);
     175           0 :     newMolecule.setForce(force);
     176           0 :     newMolecule.setPotentialEnergy(potentialEnergy);
     177             :     // add molecule to MoleculeService and LinkedCellService
     178           0 :     simplemd::Molecule* myMolecule = _moleculeService.addMolecule(newMolecule);
     179           0 :     tarch::la::Vector<MD_DIM, unsigned int> linkedCellIndex(getLinkedCellIndexForMoleculePosition(myMolecule->getConstPosition()));
     180           0 :     _linkedCellService.addMoleculeToLinkedCell(*myMolecule, linkedCellIndex);
     181           0 :   }
     182             : 
     183           0 :   tarch::la::Vector<MD_DIM, unsigned int> getLinkedCellIndexForMoleculePosition(const tarch::la::Vector<MD_DIM, double>& position) {
     184           0 :     tarch::la::Vector<MD_DIM, double> domainOffset = _parallelTopologyService.getGlobalDomainOffset();
     185           0 :     tarch::la::Vector<MD_DIM, double> meshWidth = _parallelTopologyService.getMeshWidth();
     186           0 :     tarch::la::Vector<MD_DIM, unsigned int> globalIndexOfFirstCell = _parallelTopologyService.getGlobalIndexOfFirstCell();
     187           0 :     tarch::la::Vector<MD_DIM, unsigned int> localIndexOfFirstCell = _linkedCellService.getLocalIndexOfFirstCell();
     188             : 
     189             : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
     190           0 :     const tarch::la::Vector<MD_DIM, double> domainSize = _parallelTopologyService.getGlobalDomainSize();
     191           0 :     for (unsigned int d = 0; d < MD_DIM; d++) {
     192           0 :       if ((position[d] < domainOffset[d] - meshWidth[d]) || (position[d] > domainOffset[d] + domainSize[d] + meshWidth[d])) {
     193           0 :         std::cout << "ERROR coupling::interface::impl::MDSolverInterface::addMoleculeToLinkedCell: Position ";
     194           0 :         std::cout << d << " is out of range!" << std::endl;
     195           0 :         std::cout << "Position: " << position << std::endl;
     196           0 :         exit(EXIT_FAILURE);
     197             :       }
     198             :     }
     199             : #endif
     200           0 :     tarch::la::Vector<MD_DIM, unsigned int> cellVectorIndex(0);
     201             : 
     202             :     // determine current cell index (in serial, i.e. 1-D, form)
     203           0 :     for (unsigned int d = 0; d < MD_DIM; d++) {
     204             :       // find global cell index
     205           0 :       int index = static_cast<int>(floor((position[d] - domainOffset[d]) / meshWidth[d]));
     206             :       // shift into local cell index
     207           0 :       index = index - globalIndexOfFirstCell[d] + localIndexOfFirstCell[d];
     208             : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
     209           0 :       if (index < 0) {
     210           0 :         std::cout << "ERROR coupling::interfaces::impl::SimpleMD::SimpleMDSolverInterface: index < 0: index=";
     211           0 :         std::cout << index << std::endl;
     212           0 :         exit(EXIT_FAILURE);
     213             :       }
     214             : #endif
     215           0 :       cellVectorIndex[d] = static_cast<unsigned int>(index);
     216             :     }
     217             : 
     218           0 :     return cellVectorIndex;
     219             :   }
     220             : 
     221           0 :   void setupPotentialEnergyLandscape(const tarch::la::Vector<MD_DIM, unsigned int>& indexOfFirstCouplingCell,
     222             :                                      const tarch::la::Vector<MD_DIM, unsigned int>& rangeCouplingCells,
     223             :                                      const tarch::la::Vector<MD_DIM, unsigned int>& linkedCellsPerCouplingCell) {
     224           0 :     simplemd::cellmappings::ResetPotentialEnergyMapping resetPotentialEnergyMapping;
     225           0 :     simplemd::cellmappings::LennardJonesPotentialEnergyMapping potentialEnergyMapping(_molecularPropertiesService);
     226           0 :     tarch::la::Vector<MD_DIM, unsigned int> rangeLinkedCellsExtended(0);
     227           0 :     tarch::la::Vector<MD_DIM, unsigned int> firstLinkedCell(0);
     228             :     // compute coordinates of the first linked cell and the range of linked cells to be considered
     229           0 :     for (unsigned int d = 0; d < MD_DIM; d++) {
     230             : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
     231             :       if (indexOfFirstCouplingCell[d] == 0) {
     232             :         std::cout << "ERROR setupPotentialEnergyLandscape: " << indexOfFirstCouplingCell[d] << std::endl;
     233             :         exit(EXIT_FAILURE);
     234             :       }
     235             : #endif
     236             : 
     237           0 :       firstLinkedCell[d] = (indexOfFirstCouplingCell[d] - 1) * linkedCellsPerCouplingCell[d];
     238             :       // we need to loop over one more layer of linked cells since firstLinkedCell(d) already starts one linked cell earlier
     239             :       // (e.g. already in linked cell ghost layer)
     240           0 :       rangeLinkedCellsExtended[d] = rangeCouplingCells[d] * linkedCellsPerCouplingCell[d] + _linkedCellService.getLocalIndexOfFirstCell()[d];
     241             :     }
     242             : 
     243             :     // reset potential energy first for the molecules in all relevant linked cells
     244           0 :     _linkedCellService.iterateCells(resetPotentialEnergyMapping, firstLinkedCell, rangeLinkedCellsExtended, false);
     245             : 
     246             :     // compute potential energy in these cells
     247           0 :     _linkedCellService.iterateCellPairs(potentialEnergyMapping, firstLinkedCell, rangeLinkedCellsExtended, false);
     248           0 :   }
     249             : 
     250           0 :   void calculateForceAndEnergy(coupling::interface::Molecule<MD_DIM>& molecule) {
     251           0 :     const tarch::la::Vector<MD_DIM, double> position = molecule.getPosition();
     252           0 :     const tarch::la::Vector<MD_DIM, unsigned int> linkedCellIndex = getLinkedCellIndexForMoleculePosition(position);
     253             : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
     254           0 :     for (unsigned int d = 0; d < MD_DIM; d++) {
     255           0 :       if ((linkedCellIndex[d] < 1) || (linkedCellIndex[d] > _linkedCellService.getLocalNumberOfCells()[d])) {
     256           0 :         std::cout << "ERROR coupling::interface/impl/SimpleMD/SimpleMDSolverInterface::calculateForceAndEnergy(): linkedCellIndex out of range!" << std::endl;
     257           0 :         std::cout << "LinkedCellIndex=" << linkedCellIndex << std::endl;
     258           0 :         exit(EXIT_FAILURE);
     259             :       }
     260             :     }
     261             : #endif
     262           0 :     tarch::la::Vector<MD_DIM, double> forceOld(0.0);
     263           0 :     double energy = 0.0;
     264           0 :     molecule.setPotentialEnergy(0.0);
     265             : 
     266           0 :     tarch::la::Vector<MD_DIM, unsigned int> loopIndex(0);
     267             : 
     268             : // loop over all relevant linked cells and compute force and energy contributions
     269             : #if (MD_DIM > 2)
     270           0 :     for (loopIndex[2] = linkedCellIndex[2] - 1; loopIndex[2] < linkedCellIndex[2] + 2; loopIndex[2]++) {
     271             : #endif
     272             : #if (MD_DIM > 1)
     273           0 :       for (loopIndex[1] = linkedCellIndex[1] - 1; loopIndex[1] < linkedCellIndex[1] + 2; loopIndex[1]++) {
     274             : #endif
     275           0 :         for (loopIndex[0] = linkedCellIndex[0] - 1; loopIndex[0] < linkedCellIndex[0] + 2; loopIndex[0]++) {
     276             : 
     277             :           // loop over all molecules in each cell
     278           0 :           const simplemd::LinkedCell& thisCell = _linkedCellService.getLinkedCell(loopIndex);
     279           0 :           for (std::list<simplemd::Molecule*>::const_iterator it = thisCell.constBegin(); it != thisCell.constEnd(); it++) {
     280           0 :             forceOld += getLennardJonesForce(position, (*it)->getConstPosition());
     281           0 :             energy += getPotentialEnergy(position, (*it)->getConstPosition());
     282             :           }
     283             :         }
     284             : #if (MD_DIM > 1)
     285             :       }
     286             : #endif
     287             : #if (MD_DIM > 2)
     288             :     }
     289             : #endif
     290             :     // set force in forceOld-entry
     291           0 :     molecule.setForce(forceOld);
     292             :     // set total potential energy
     293           0 :     molecule.setPotentialEnergy(energy);
     294           0 :   }
     295             : 
     296           0 :   void synchronizeMoleculesAfterMassModification() {
     297           0 :     _boundaryTreatment.emptyGhostBoundaryCells();
     298             : #if (TEST_TCHIPEV == MD_NO)
     299           0 :     _boundaryTreatment.fillBoundaryCells(_localBoundaryInformation, _parallelTopologyService);
     300             : #else
     301             :     _boundaryTreatment.putBoundaryParticlesToInnerCellsAndFillBoundaryCells(_localBoundaryInformation, _parallelTopologyService);
     302             : #endif
     303           0 :   }
     304             : 
     305           0 :   void synchronizeMoleculesAfterMomentumModification() {}
     306             : 
     307           0 :   double getDt() { return _dt; }
     308             : 
     309           0 :   coupling::interface::MoleculeIterator<simplemd::LinkedCell, MD_DIM>* getMoleculeIterator(simplemd::LinkedCell& cell) {
     310           0 :     return new coupling::interface::SimpleMDMoleculeIterator(cell);
     311             :   }
     312             : };
     313             : 
     314             : } // namespace interface
     315             : } // namespace coupling
     316             : #endif // _MOLECULARDYNAMICS_COUPLING_INTERFACE_SIMPLEMDSOLVERINTERFACE_CPP_