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 16 : 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
|
