ls1mardyn / ls1-mardyn

ls1-MarDyn is a massively parallel Molecular Dynamics (MD) code for large systems. Its main target is the simulation of thermodynamics and nanofluidics. ls1-MarDyn is designed with a focus on performance and easy extensibility.
http://www.ls1-mardyn.de
Other
28 stars 15 forks source link

Fix of mixing rules #333

Closed HomesGH closed 2 months ago

HomesGH commented 2 months ago

Description

The mixing rules are kind of messed up for now, see issues:

73

74

75

331

There was already an attempt #76 for a fix but it was not merged but closed.

~It would be elegant to use a lambda function for the application of the (here and here) rule to make it more generic. Unfortunately, I wasn't able to make it work due to the scopes here. Help is appreciated 😉 The open points above would also be solved by applying the LB-specific lambda function in the Comp2Param.~ Done

FG-TUM commented 2 months ago

There was already an attempt #76 for a fix but it was not merged but closed.

It would be elegant to use a lambda function for the application of the (here and here) rule to make it more generic. Unfortunately, I wasn't able to make it work due to the scopes here. Help is appreciated 😉

Sorry I don't fully understand what is the problem here. Can you clarify?

HomesGH commented 2 months ago

There was already an attempt #76 for a fix but it was not merged but closed. It would be elegant to use a lambda function for the application of the (here and here) rule to make it more generic. Unfortunately, I wasn't able to make it work due to the scopes here. Help is appreciated 😉

Sorry I don't fully understand what is the problem here. Can you clarify?

For now, the LB rule is used in the general calculations in Comp2Params: https://github.com/ls1mardyn/ls1-mardyn/blob/4f5b49f4e42b30245fd865baceedc20099f3c1b3/src/molecules/Comp2Param.cpp#L65-L66

Since other rules might be introduced later, this is not ideal and it would be better to do something like (note the lambda expressions):

if (mixingrule->getType() == "LB") {
                eta = mixingrule->getParameters().at(0);
                xi  = mixingrule->getParameters().at(1);
                auto mixingSigma = [&](double epsi, double epsj) { return xi * sqrt(epsi * epsj); }
                auto mixingEpsilon = [&](double epsi, double epsj) { return xi * sqrt(epsi * epsj); }
// #ifndef NDEBUG
                Log::global_log->info() << "Mixing : cid+1(compi)=" << compi+1 << " <--> cid+1(compj)=" << compj+1 << ": xi=" << xi << ", eta=" << eta << std::endl;
// #endif
            } else {
                Log::global_log->error() << "Mixing: Only LB rule supported" << std::endl;
                Simulation::exit(1);
            }

and later in the calculations

epsilon24 = 24. * mixingEpsilon (epsi, epsj);
sigma2 = .5 * mixingSigma (sigi, sigj);

instead of the current (LB-specific) implementation.

But with this, the lambda expressions are just defined in the scope of the if statement.

FG-TUM commented 2 months ago

Ah ok so you want to store functions in the variables and later use them to get the actual mixing coefficients, correct?

You could do something like this:

std::function<double(double, double)> mixingEpsilon;
std::function<double(double, double)> mixingSigma;

if (mixingrule->getType() == "LB") {
    // ...
    mixingSigma = [=](double epsi, double epsj) { return xi * sqrt(epsi * epsj); };
    mixingEpsilon = [=](double epsi, double epsj) { return xi * sqrt(epsi * epsj); };
}

Note the use of taking xi and epsi by copy ([=]) since the here-referenced variables probably won't exist anymore when you want to invoke the lambdas.

I've almost never used the std::function interface. Therefore, I have no idea about its performance implications. Where exactly would you want to invoke the functions?

HomesGH commented 2 months ago

I've almost never used the std::function interface. Therefore, I have no idea about its performance implications. Where exactly would you want to invoke the functions?

Thank you! I have used std::function now. The code is used during the initialization of the comp2params and, therefore, just called once or twice at the beginning of the simulation.