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_CELLMAPPINGS_ZHOUBOUNDARYFORCE_H_
6 : #define _MOLECULARDYNAMICS_COUPLING_CELLMAPPINGS_ZHOUBOUNDARYFORCE_H_
7 :
8 : #include "coupling/interface/MDSolverInterface.h"
9 : #include "coupling/interface/Molecule.h"
10 : #include "coupling/interface/MoleculeIterator.h"
11 : #include <cmath>
12 : #include <iostream>
13 :
14 : namespace coupling {
15 : namespace cellmappings {
16 : template <class LinkedCell, unsigned int dim> class ZhouBoundaryForce;
17 : }
18 : } // namespace coupling
19 :
20 : /** applies the Zhou boundary force to all molecules assuming a cut-off radius
21 : *r_c=2.5. We only consider molecules in the outermost coupling cell. Thus,
22 : *the current method will only be employed to molecules within a
23 : *distance=min(r_c,size of coupling cell). The force is based on the research
24 : *paper: W.J. Zhou, H.B. Luan, Y.L. He, J. Sun, W.Q. Tao. A study on boundary
25 : *force model used in multiscale simulations with non-periodic boundary
26 : *condition Microfluid Nanofluid 16(3): 587-595, 2014
27 : *
28 : * The force is computed based on an interpolation formula which has been
29 : *obtained from various MD simulation runs at different temperature and density
30 : *values. The investigated range is:
31 : * 0.4*m*sigma-3<density<0.9*m*sigma-3, 1.3*eps/kB<temperature<3.9*eps/kB.
32 : * Moreover, only 3D simulations have been considered in this study.
33 : * @brief This class applies the Zhou boundary force to all molecules
34 : *assuming a cut-off radius r_c=2.5.
35 : * @tparam LinkedCell cell type
36 : * @tparam dim Number of dimensions; it can be 1, 2 or 3
37 : * @author Philipp Neumann
38 : * \todo Philipp please take a look on this class
39 : */
40 : template <class LinkedCell, unsigned int dim> class coupling::cellmappings::ZhouBoundaryForce {
41 : public:
42 : /** currently, we support constant temperature and density and 3D only. To be
43 : *extended
44 : * @param density
45 : * @param temperature
46 : * @param epsilon
47 : * @param sigma
48 : * @param boundary
49 : * @param domainOffset
50 : * @param domainSize
51 : * @param mdSolverInterface
52 : */
53 0 : ZhouBoundaryForce(const double& density, const double& temperature, const double& epsilon, const double& sigma,
54 : const tarch::la::Vector<2 * dim, bool>& boundary, const tarch::la::Vector<dim, double>& domainOffset,
55 : const tarch::la::Vector<dim, double>& domainSize, coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface)
56 0 : : _boundary(boundary), _domainLowerLeft(domainOffset), _domainUpperRight(domainSize + domainOffset), _mdSolverInterface(mdSolverInterface),
57 0 : _p1(getP1(density, temperature)), _p2(getP2(density, temperature)), _p3(getP3(density, temperature)), _q1(getQ1(density)), _q2(getQ2(density)),
58 0 : _q3(getQ3(density)), _forceFactor(epsilon / sigma), _sigma(sigma),
59 : /** Fixed step width for numerical integration of boundary force term */
60 0 : _energyResolution(1000), _energyTable(new double[_energyResolution]) {
61 : // simple trapezoidal integration scheme
62 0 : double force(0);
63 0 : double energy(0);
64 0 : double dx(2.5 * sigma / (_energyResolution - 1));
65 0 : for (unsigned int i = 0; i < _energyResolution; i++) {
66 0 : force = getScalarBoundaryForce(2.5 - i * dx);
67 0 : _energyTable[i] = energy + force / 2 * dx;
68 0 : energy += force * dx;
69 : }
70 :
71 : // for(unsigned int i=0;i<_energyResolution;i++){
72 : // std::cout << "_energyTable[" << i << "] = " << _energyTable[i] <<
73 : // std::endl;
74 : // }
75 : // for(double pos = 0; pos <= 2.5; pos += .25/_energyResolution){
76 : // std::cout << "getPotentialEnergy(" << pos << ") = " <<
77 : // getPotentialEnergy(pos) << std::endl;
78 : // }
79 0 : }
80 :
81 : /** Destructor */
82 0 : ~ZhouBoundaryForce() { delete[] _energyTable; }
83 :
84 : /** empty function
85 : */
86 : void beginCellIteration() {}
87 :
88 : /** empty function
89 : */
90 : void endCellIteration() {}
91 :
92 : /** extracts position and force of each molecule, add boundary force and
93 : *applies it to Molecule
94 : * @param cell
95 : */
96 0 : void handleCell(LinkedCell& cell) {
97 0 : coupling::interface::MoleculeIterator<LinkedCell, dim>* it = _mdSolverInterface->getMoleculeIterator(cell);
98 0 : it->begin();
99 0 : while (it->continueIteration()) {
100 0 : coupling::interface::Molecule<dim>& wrapper(it->get());
101 : // std::cout << "apply force to Molecule " << wrapper.getPosition() <<
102 : // std::endl;
103 :
104 : // extract position and force of each molecule
105 0 : const tarch::la::Vector<dim, double> position(wrapper.getPosition());
106 0 : tarch::la::Vector<dim, double> force(wrapper.getForce());
107 : // add boundary force
108 0 : force = force + getBoundaryForces(position);
109 0 : wrapper.setForce(force);
110 :
111 0 : it->next();
112 : }
113 0 : delete it;
114 0 : }
115 :
116 : /** returns the potential energy
117 : * @param position
118 : * @return energy
119 : */
120 0 : double getPotentialEnergy(const tarch::la::Vector<dim, double>& position) const {
121 0 : double energy(0);
122 0 : const double distance = 2.5 * _sigma;
123 0 : for (unsigned int d = 0; d < dim; d++) {
124 0 : if (_boundary[2 * d] && (position[d] - _domainLowerLeft[d] < distance)) {
125 0 : energy += getPotentialEnergy((position[d] - _domainLowerLeft[d]) / _sigma);
126 : }
127 0 : if (_boundary[2 * d + 1] && (_domainUpperRight[d] - position[d] < distance)) {
128 0 : energy += getPotentialEnergy((_domainUpperRight[d] - position[d]) / _sigma);
129 : }
130 : }
131 0 : return energy;
132 : }
133 :
134 : /** checks the boundary flags for open boundaries. If an open boundary is
135 : *encountered and this molecule is close to that boundary (distance smaller
136 : *than 2.5, cf. Zhou paper), the respective force contribution in that
137 : *dimension is added/subtracted.
138 : * @param position
139 : * @return force
140 : */
141 0 : tarch::la::Vector<dim, double> getBoundaryForces(const tarch::la::Vector<dim, double>& position) const {
142 0 : tarch::la::Vector<dim, double> force(0.0);
143 0 : const double distance = 2.5 * _sigma;
144 :
145 0 : for (unsigned int d = 0; d < dim; d++) {
146 : // for left/lower/back boundary: add force; for right/upper/front
147 : // boundary: subtract force value
148 0 : if (_boundary[2 * d] && (position[d] - _domainLowerLeft[d] < distance)) {
149 0 : force[d] += _forceFactor * getScalarBoundaryForce((position[d] - _domainLowerLeft[d]) / _sigma);
150 : }
151 0 : if (_boundary[2 * d + 1] && (_domainUpperRight[d] - position[d] < distance)) {
152 0 : force[d] -= _forceFactor * getScalarBoundaryForce((_domainUpperRight[d] - position[d]) / _sigma);
153 : }
154 : }
155 0 : return force;
156 : }
157 :
158 : private:
159 : /** perform interpolation between nearest energy lookup table values
160 : * @param rw
161 : * @return nearest energy to the ????
162 : */
163 0 : double getPotentialEnergy(double rw) const {
164 0 : double dx(2.5 * _sigma / (_energyResolution - 1));
165 0 : int upper = (int)((2.5 - rw) / 2.5 * (_energyResolution - 1) - .5);
166 0 : int lower = (int)((2.5 - rw) / 2.5 * (_energyResolution - 1) + .5);
167 0 : return _energyTable[lower] * ((2.5 - upper * dx - rw) / dx) + _energyTable[upper] * ((rw - 2.5 + lower * dx) / dx);
168 : }
169 :
170 : /** evaluates the force expression for a scalar component of the force vector.
171 : *We thus add up dimensional contributions if several boundaries are located
172 : *beside each other.
173 : * @param rw
174 : * @return Scalar boundary force
175 : */
176 0 : double getScalarBoundaryForce(double rw) const {
177 0 : if (rw < 1.04) {
178 0 : return _p1 + _p2 * exp(pow((rw + 0.25), 3.4)) * cos(_p3 * rw);
179 : } else {
180 0 : return -1.0 / (_q1 + _q2 * (2.5 - rw) * (2.5 - rw) + _q3 / ((2.5 - rw) * (2.5 - rw)));
181 : }
182 : }
183 :
184 : /** computes and returns the coefficient p1 according to Zhou-paper
185 : * @param density
186 : * @param temperature
187 : * @return p1
188 : */
189 0 : double getP1(const double& density, const double& temperature) const {
190 0 : const double logrho = log(density);
191 0 : const double T2 = temperature * temperature;
192 0 : const double T3 = T2 * temperature;
193 :
194 : // an additional factor 10**-5 is in the expression 4.599/(...) -> private
195 : // communication with Wenjing Zhou (typo in paper)
196 0 : const double p1 = (-18.953 + 53.369 * temperature - 1.253 * T2 + 4.599 / 100000.0 * T3 + 59.871 * logrho + 19.737 * logrho * logrho) /
197 0 : (1.0 + 2.592 * temperature - 0.557 * T2 + 0.049 * T3 - 13.912 * logrho + 18.657 * logrho * logrho);
198 0 : return p1;
199 : }
200 :
201 : /** computes and returns the coefficient p2 according to Zhou-paper
202 : * @param density
203 : * @param temperature
204 : * @return p2
205 : */
206 0 : double getP2(const double& density, const double& temperature) const {
207 0 : const double logrho = log(density);
208 0 : const double T2 = temperature * temperature;
209 0 : const double T3 = T2 * temperature;
210 :
211 0 : const double p2 = (-0.094 + 2.808 * temperature - 0.019 * T2 - 0.001 * T3 + 2.823 * logrho + 2.071 * logrho * logrho) /
212 0 : (1.0 + 0.168 * temperature - 0.013 * T2 - 4.323 * logrho + 2.557 * logrho * logrho - 2.155 * logrho * logrho * logrho);
213 0 : return p2;
214 : }
215 :
216 : /** computes and returns the coefficient p3 according to Zhou-paper
217 : * @param density
218 : * @param temperature
219 : * @return p3
220 : */
221 0 : double getP3(const double& density, const double& temperature) const {
222 0 : const double T394 = pow(temperature, 0.394);
223 0 : const double rho17437 = pow(density, 17.437);
224 0 : const double p3 = 3.934 + 0.099 * T394 - 0.097 * rho17437 + 0.075 * T394 * rho17437;
225 0 : return p3;
226 : }
227 :
228 : /** computes and returns the coefficient q1 according to Zhou-paper
229 : * @param density
230 : * @return q1
231 : */
232 0 : double getQ1(const double& density) const {
233 0 : const double rho2 = density * density;
234 0 : const double rho3 = rho2 * density;
235 0 : const double rho4 = rho2 * rho2;
236 0 : return -30.471 + 113.879 * density - 207.205 * rho2 + 184.242 * rho3 - 62.879 * rho4;
237 : }
238 :
239 : /** computes and returns the coefficient q2 according to Zhou-paper
240 : * @param density
241 : * @return q2
242 : */
243 0 : double getQ2(const double& density) const {
244 0 : const double rho2 = density * density;
245 0 : const double rho3 = rho2 * density;
246 0 : const double rho4 = rho2 * rho2;
247 0 : return 6.938 - 25.788 * density + 46.773 * rho2 - 41.768 * rho3 + 14.394 * rho4;
248 : }
249 :
250 : /** computes and returns the coefficient q3 according to Zhou-paper
251 : * @param density
252 : * @return q3
253 : */
254 0 : double getQ3(const double& density) const {
255 0 : const double rho2 = density * density;
256 0 : const double rho3 = rho2 * density;
257 0 : const double rho4 = rho2 * rho2;
258 0 : return 39.634 - 147.821 * density + 269.519 * rho2 - 239.066 * rho3 + 81.439 * rho4;
259 : }
260 :
261 : /** an entry is true, if this is an open boundary. For enumeration of
262 : * boundaries, see BoundaryForceConfiguration */
263 : const tarch::la::Vector<2 * dim, bool> _boundary;
264 : /** lower left corner of MD domain */
265 : const tarch::la::Vector<dim, double> _domainLowerLeft;
266 : /** upper right corner of MD domain */
267 : const tarch::la::Vector<dim, double> _domainUpperRight;
268 :
269 : coupling::interface::MDSolverInterface<LinkedCell, dim>* const _mdSolverInterface;
270 : // helper variables to speed up evaluation, see Zhou paper
271 : const double _p1;
272 : const double _p2;
273 : const double _p3;
274 : const double _q1;
275 : const double _q2;
276 : const double _q3;
277 :
278 : /** factor to scale the force according to MD units */
279 : const double _forceFactor;
280 : /** factor to scale length according to MD units */
281 : const double _sigma;
282 :
283 : const unsigned int _energyResolution;
284 : double* _energyTable;
285 : };
286 : #endif // _MOLECULARDYNAMICS_COUPLING_CELLMAPPINGS_ZHOUBOUNDARYFORCE_H_
|
