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             : 
       6             : namespace coupling {
       7             : namespace interface {
       8             : /** defines a MDsimulation based on simplemd::MolecularDynamicsSimulation but for coupled MD simulations.
       9             :  * Thus, the implementation of one timestep slightly differs from the one of the base class.
      10             :  *      @brief defines MD simulation from default simple MD code.
      11             :  *  @author Philipp Neumann
      12             :  */
      13             : class SimpleMDSimulation : public coupling::interface::MDSimulation, public simplemd::MolecularDynamicsSimulation {
      14             : public:
      15          12 :   SimpleMDSimulation(const simplemd::configurations::MolecularDynamicsConfiguration& configuration)
      16          12 :       : coupling::interface::MDSimulation(), simplemd::MolecularDynamicsSimulation(configuration), _couplingCellService(NULL), _couplingSwitchedOn(true) {}
      17             : 
      18          12 :   virtual ~SimpleMDSimulation() {}
      19             : 
      20          16 :   virtual void switchOffCoupling() { _couplingSwitchedOn = false; }
      21             : 
      22           8 :   virtual void switchOnCoupling() { _couplingSwitchedOn = true; }
      23             : 
      24          24 :   virtual void simulateTimesteps(const unsigned int& numberTimesteps, const unsigned int& firstTimestep) {
      25         144 :     for (unsigned int t = firstTimestep; t < firstTimestep + numberTimesteps; t++) {
      26         120 :       simulateOneCouplingTimestep(t);
      27             :     }
      28          24 :   }
      29             : 
      30         120 :   void simulateOneCouplingTimestep(const unsigned int& t) {
      31             :     // if coupling is switched off, perform "normal" MD timestep
      32         120 :     if (!_couplingSwitchedOn) {
      33         120 :       simulateOneTimestep(t);
      34         120 :       return;
      35             :     }
      36           0 :     if (_parallelTopologyService->isIdle()) {
      37             :       return;
      38             :     }
      39             : 
      40           0 :     _boundaryTreatment->putBoundaryParticlesToInnerCellsAndFillBoundaryCells(_localBoundary, *_parallelTopologyService);
      41             : 
      42             :     // call to synchronise data in cells; needs to be at this point of the
      43             :     // coupling algorithm as the particles need to be placed inside the correct
      44             :     // sampling volumes (hence: after communication with neighbours and molecule
      45             :     // updates); do it BEFORE quantities are manipulated as we can then also do
      46             :     // some pre-processing here.
      47           0 :     _couplingCellService->processInnerCouplingCellAfterMDTimestep();
      48             :     // ------------ coupling step: distribute mass ---------------------
      49           0 :     _couplingCellService->distributeMass(t);
      50             :     // for isothermal simulations: apply thermostat
      51           0 :     _couplingCellService->applyTemperatureToMolecules(t);
      52             : 
      53             :     // ---------- from here: go on with usual MD algorithm
      54             :     // ------------------------------
      55             : 
      56             :     // compute forces. After this step, each molecule has received all force
      57             :     // contributions from its neighbors.
      58           0 :     _linkedCellService->iterateCellPairs(*_lennardJonesForce, false);
      59             : 
      60             :     // distribute momentum -> some methods require modification of force terms,
      61             :     // therefore we call it AFTER the force computation and before everything else
      62           0 :     _couplingCellService->distributeMomentum(t);
      63             :     // apply boundary forces
      64           0 :     _couplingCellService->applyBoundaryForce(t);
      65             :     // evaluate statistics
      66           0 :     evaluateStatistics(t);
      67             : 
      68           0 :     _boundaryTreatment->emptyGhostBoundaryCells();
      69             : 
      70             :     // plot VTK output
      71           0 :     if ((_configuration.getVTKConfiguration().getWriteEveryTimestep() != 0) && (t % _configuration.getVTKConfiguration().getWriteEveryTimestep() == 0)) {
      72           0 :       _vtkMoleculeWriter->setTimestep(t);
      73           0 :       _moleculeService->iterateMolecules(*_vtkMoleculeWriter, false);
      74             :     }
      75             : 
      76             : // plot ADIOS2 output
      77             : #if BUILD_WITH_ADIOS2
      78             :     if ((_configuration.getAdios2Configuration().getWriteEveryTimestep() != 0) && (t % _configuration.getAdios2Configuration().getWriteEveryTimestep() == 0)) {
      79             :       _Adios2Writer->setTimestep(t);
      80             :       _moleculeService->iterateMolecules(*_Adios2Writer, false);
      81             :     }
      82             : #endif
      83             : 
      84             :     // write checkpoint
      85           0 :     if ((_configuration.getCheckpointConfiguration().getWriteEveryTimestep() != 0) &&
      86           0 :         (t % _configuration.getCheckpointConfiguration().getWriteEveryTimestep() == 0)) {
      87           0 :       _moleculeService->writeCheckPoint(*_parallelTopologyService, _configuration.getCheckpointConfiguration().getFilename(), t);
      88             :     }
      89             :     // reorganise memory if needed
      90           0 :     if ((_configuration.getSimulationConfiguration().getReorganiseMemoryEveryTimestep() != 0) &&
      91           0 :         (t % _configuration.getSimulationConfiguration().getReorganiseMemoryEveryTimestep() == 0)) {
      92           0 :       _moleculeService->reorganiseMemory(*_parallelTopologyService, *_linkedCellService);
      93             :     }
      94             :     // plot also coupling cell information
      95           0 :     _couplingCellService->plotEveryMicroscopicTimestep(t);
      96             : 
      97           0 :     _linkedCellService->iterateCells(*_emptyLinkedListsMapping, false);
      98             : 
      99             :     // time integration. After this step, the velocities and the positions of the
     100             :     // molecules have been updated.
     101           0 :     _moleculeService->iterateMolecules(*_timeIntegrator, false);
     102             : 
     103             :     // sort molecules into linked cells
     104           0 :     _moleculeService->iterateMolecules(*_updateLinkedCellListsMapping, false);
     105             : 
     106           0 :     if (_parallelTopologyService->getProcessCoordinates() == tarch::la::Vector<MD_DIM, unsigned int>(0)) {
     107             :       // if(t%50==0) std::cout <<"Finish MD timestep " << t << "..." << std::endl;
     108             :     }
     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 :     getMoleculeService().writeCheckPoint(getParallelTopologyService(), filestem, t);
     134           4 :   }
     135             : 
     136             :   // function particularly needed to init MD solver interface -> should only be
     137             :   // called from factory
     138          12 :   simplemd::BoundaryTreatment& getBoundaryTreatment() { return *_boundaryTreatment; }
     139          28 :   simplemd::services::ParallelTopologyService& getParallelTopologyService() { return *_parallelTopologyService; }
     140          16 :   simplemd::services::MoleculeService& getMoleculeService() {
     141             :     #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
     142          16 :     if(_moleculeService == NULL){
     143           0 :       std::cout <<"ERROR coupling::interface::MDSimulation::getMoleculeService(): _moleculeService == NULL " << std::endl;
     144           0 :       exit(1);
     145             :     }
     146             :     #endif
     147          16 :     return *_moleculeService;
     148             :   }
     149          12 :   simplemd::services::LinkedCellService& getLinkedCellService() { return *_linkedCellService; }
     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             : };
     167             : } // namespace interface
     168             : } // namespace coupling