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