Line data    Source code

       1             : #pragma once
       2             : #include "coupling/interface/MDSimulation.h"
       3             : #include "coupling/interface/MamicoInterfaceProvider.h"
       4             : #include "simplemd/MolecularDynamicsSimulation.h"
       5             : #include "simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h"
       6             : #include "simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h"
       7             : #include "simplemd/molecule-mappings/WriteCheckPointMapping.h"
       8             : 
       9             : namespace coupling {
      10             : namespace interface {
      11             : /** defines a MDsimulation based on simplemd::MolecularDynamicsSimulation but for coupled MD simulations.
      12             :  * Thus, the implementation of one timestep slightly differs from the one of the base class.
      13             :  *      @brief defines MD simulation from default simple MD code.
      14             :  *  @author Philipp Neumann
      15             :  */
      16             : class SimpleMDSimulation : public coupling::interface::MDSimulation, public simplemd::MolecularDynamicsSimulation {
      17             : public:
      18          12 :   SimpleMDSimulation(const simplemd::configurations::MolecularDynamicsConfiguration& configuration)
      19          12 :       : coupling::interface::MDSimulation(), simplemd::MolecularDynamicsSimulation(configuration), _couplingCellService(NULL), _couplingSwitchedOn(true) {}
      20             : 
      21          12 :   virtual ~SimpleMDSimulation() {}
      22             : 
      23          16 :   virtual void switchOffCoupling() { _couplingSwitchedOn = false; }
      24             : 
      25           8 :   virtual void switchOnCoupling() { _couplingSwitchedOn = true; }
      26             : 
      27          24 :   virtual void simulateTimesteps(const unsigned int& numberTimesteps, const unsigned int& firstTimestep) {
      28         144 :     for (unsigned int t = firstTimestep; t < firstTimestep + numberTimesteps; t++) {
      29         120 :       simulateOneCouplingTimestep(t);
      30             :     }
      31          24 :   }
      32             : 
      33         120 :   void simulateOneCouplingTimestep(const unsigned int& t) {
      34             :     // if coupling is switched off, perform "normal" MD timestep
      35         120 :     if (!_couplingSwitchedOn) {
      36         120 :       simulateOneTimestep(t);
      37         120 :       return;
      38             :     }
      39           0 :     if (_parallelTopologyService->isIdle()) {
      40             :       return;
      41             :     }
      42             : 
      43           0 :     _boundaryTreatment->putBoundaryParticlesToInnerCellsAndFillBoundaryCells(_localBoundary, *_parallelTopologyService);
      44             : 
      45             :     // call to synchronise data in cells; needs to be at this point of the
      46             :     // coupling algorithm as the particles need to be placed inside the correct
      47             :     // sampling volumes (hence: after communication with neighbours and molecule
      48             :     // updates); do it BEFORE quantities are manipulated as we can then also do
      49             :     // some pre-processing here.
      50           0 :     _couplingCellService->processInnerCouplingCellAfterMDTimestep();
      51             :     // ------------ coupling step: distribute mass ---------------------
      52           0 :     _couplingCellService->distributeMass(t);
      53             :     // for isothermal simulations: apply thermostat
      54           0 :     _couplingCellService->applyTemperatureToMolecules(t);
      55             : 
      56             :     // ---------- from here: go on with usual MD algorithm
      57             :     // ------------------------------
      58             : 
      59             :     // compute forces. After this step, each molecule has received all force
      60             :     // contributions from its neighbors.
      61             : #if (TARCH_DEBUG == TARCH_YES)
      62           0 :     _moleculeService->getContainer().iterateMolecules(_convertForcesFloatToFixedMapping);
      63             : #endif
      64           0 :     _moleculeService->getContainer().iterateCellPairs(*_lennardJonesForce);
      65             : #if (TARCH_DEBUG == TARCH_YES)
      66           0 :     _moleculeService->getContainer().iterateMolecules(_convertForcesFixedToFloatMapping);
      67             : #endif
      68             : 
      69             :     // distribute momentum -> some methods require modification of force terms,
      70             :     // therefore we call it AFTER the force computation and before everything else
      71           0 :     _couplingCellService->distributeMomentum(t);
      72             :     // apply boundary forces
      73           0 :     _couplingCellService->applyBoundaryForce(t);
      74             :     // evaluate statistics
      75           0 :     evaluateStatistics(t);
      76             : 
      77           0 :     _boundaryTreatment->emptyGhostBoundaryCells();
      78             : 
      79             :     // plot VTK output
      80           0 :     if ((_configuration.getVTKConfiguration().getWriteEveryTimestep() != 0) && (t % _configuration.getVTKConfiguration().getWriteEveryTimestep() == 0)) {
      81           0 :       _vtkMoleculeWriter->setTimestep(t);
      82           0 :       _moleculeService->getContainer().iterateMolecules(*_vtkMoleculeWriter);
      83             :     }
      84             : 
      85             : // plot ADIOS2 output
      86             : #if BUILD_WITH_ADIOS2
      87             :     if ((_configuration.getAdios2Configuration().getWriteEveryTimestep() != 0) && (t % _configuration.getAdios2Configuration().getWriteEveryTimestep() == 0)) {
      88             :       _Adios2Writer->setTimestep(t);
      89             :       _moleculeService->getContainer().iterateMolecules(*_Adios2Writer);
      90             :     }
      91             : #endif
      92             : 
      93             :     // write checkpoint
      94           0 :     if ((_configuration.getCheckpointConfiguration().getWriteEveryTimestep() != 0) &&
      95           0 :         (t % _configuration.getCheckpointConfiguration().getWriteEveryTimestep() == 0)) {
      96           0 :       simplemd::moleculemappings::WriteCheckPointMapping writeCheckPointMapping(*_parallelTopologyService,
      97           0 :                                                                                 _configuration.getCheckpointConfiguration().getFilename(), t);
      98           0 :       _moleculeService->getContainer().iterateMolecules(writeCheckPointMapping);
      99           0 :     }
     100             :     // plot also coupling cell information
     101           0 :     _couplingCellService->plotEveryMicroscopicTimestep(t);
     102             : 
     103             :     // time integration. After this step, the velocities and the positions of the
     104             :     // molecules have been updated.
     105           0 :     _moleculeService->getContainer().iterateMolecules(*_timeIntegrator);
     106             : 
     107             :     // sort molecules into linked cells
     108           0 :     _moleculeService->getContainer().sort();
     109             :   }
     110             : 
     111           0 :   virtual void sortMoleculesIntoCells() {
     112             :     // nop required, since the linked cells are very tightly linked to mamico
     113           0 :   }
     114             : 
     115           0 :   virtual void setCouplingCellService(coupling::services::CouplingCellService<MDSIMULATIONFACTORY_DIMENSION>* couplingCellService) {
     116           0 :     _couplingCellService = couplingCellService;
     117             :     // set the cell service also in singleton of mamico interface provider ->
     118             :     // typically not required in coupling, but makes the simulation state more
     119             :     // consistent compared to using LAMMPS
     120           0 :     coupling::interface::MamicoInterfaceProvider<simplemd::LinkedCell, MDSIMULATIONFACTORY_DIMENSION>::getInstance().setCouplingCellService(
     121             :         couplingCellService);
     122           0 :   }
     123             : 
     124           0 :   virtual void init() { initServices(); }
     125             : 
     126          12 :   virtual void init(const tarch::utils::MultiMDService<MDSIMULATIONFACTORY_DIMENSION>& multiMDService, unsigned int localMDSimulation) {
     127          12 :     initServices(multiMDService, localMDSimulation);
     128          12 :   }
     129             : 
     130          12 :   virtual void shutdown() { shutdownServices(); }
     131             : 
     132           4 :   virtual void writeCheckpoint(const std::string& filestem, const unsigned int& t) {
     133           4 :     simplemd::moleculemappings::WriteCheckPointMapping writeCheckPointMapping(getParallelTopologyService(), filestem, t);
     134           4 :     _moleculeService->getContainer().iterateMolecules(writeCheckPointMapping);
     135           4 :   }
     136             : 
     137             :   // function particularly needed to init MD solver interface -> should only be
     138             :   // called from factory
     139          12 :   simplemd::BoundaryTreatment& getBoundaryTreatment() { return *_boundaryTreatment; }
     140          28 :   simplemd::services::ParallelTopologyService& getParallelTopologyService() { return *_parallelTopologyService; }
     141          12 :   simplemd::services::MoleculeService& getMoleculeService() {
     142             : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
     143          12 :     if (_moleculeService == NULL) {
     144           0 :       std::cout << "ERROR coupling::interface::MDSimulation::getMoleculeService(): _moleculeService == NULL " << std::endl;
     145           0 :       exit(1);
     146             :     }
     147             : #endif
     148          12 :     return *_moleculeService;
     149             :   }
     150          12 :   const simplemd::services::MolecularPropertiesService& getMolecularPropertiesService() {
     151             : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
     152          12 :     if (_molecularPropertiesService == NULL) {
     153           0 :       std::cout << "ERROR coupling::interface::MDSimulation::getMolecularPropertiesService(): _molecularPropertiesService == NULL " << std::endl;
     154           0 :       exit(1);
     155             :     }
     156             : #endif
     157          12 :     return *_molecularPropertiesService;
     158             :   }
     159             : 
     160             : private:
     161             :   /** @brief the coupling cell service for the coupled md simulation  */
     162             :   coupling::services::CouplingCellService<MD_DIM>* _couplingCellService;
     163             :   /** @brief bool holding the current state of the coupling: true - coupled
     164             :    * simulation and false - independent md simulation */
     165             :   bool _couplingSwitchedOn;
     166             : #if (TARCH_DEBUG == TARCH_YES)
     167             :   simplemd::moleculemappings::ConvertForcesFloatToFixedMapping _convertForcesFloatToFixedMapping;
     168             :   simplemd::moleculemappings::ConvertForcesFixedToFloatMapping _convertForcesFixedToFloatMapping;
     169             : #endif
     170             : };
     171             : } // namespace interface
     172             : } // namespace coupling