gianni/hybridmd - Multiscale Laboratory

Hydrodynamics-molecular coupling

This is a work in progress page to expose the main components of the coupling molecular hydrodynamics code of references:

R. Delgado-Buscalioni and G. De Fabritiis, Embedding molecular dynamics within fluctuating hydrodynamics in multiscale simulations of liquids, Phys. Rev. E 76, 036709 (2007).
G. De Fabritiis, M. Serrano, R. Delgado-Buscalioni and P. V. Coveney, Fluctuating hydrodynamic modelling of fluids at the nanoscale, Phys. Rev. E 75, 026307 (2007).
G. De Fabritiis, R. Delgado-Buscalioni and P. V. Coveney, Modelling the mesoscale with molecular specificity, Phys. Rev. Lett. 97, 134501 (2006).
G. Giupponi, G. De Fabritiis and P. V. Coveney, Hybrid method coupling fluctuating hydrodynamics and molecular dynamics for the simulation of macromolecules , J. Chem. Phys. 126, 154903 (2007).

All the part codes in this page are under the GPL3 license, which you must accept in order to use the code.

Main coupling routine to apply a Stokes drag coming from an hydrodynamic code into the molecular system and vice-versa

void interpolator(Mesofluid& meso, const vector3& pos, double weights[8], int cells[8]) {
  int x,y,z;
  vector3 length = meso.Delta;

  meso.find_box(pos - 0.5*length , x, y, z);//center in dual space
  vector3 O = meso.center_point(x, y, z);
  vector3 relpos = meso.delta(pos, O);


  relpos.x() /= length.x();
  relpos.y() /= length.y();
  relpos.z() /= length.z();

  double W=0.0;
  for(int i = 0; i <= 1; i++)
  for(int j = 0; j <= 1; j++)
  for(int k = 0; k <= 1; k++) {
    int whoami = 4*i + 2*j + k;
    double weight =  fabs((1 - i - relpos.x())
            * (1 - j - relpos.y())
            * (1 - k - relpos.z()));
    weights[ whoami ] = weight;
    W += weight;
    cells[ whoami ] = meso.index(x + i, y + j, z + k);
  }
  if (fabs(W-1)>10e-4) cout<<"Debug: W-1 "<<W-1<<" "<<relpos<<" "<<pos<<" x y z "<<x<<' '<<y<<' '<<z<<endl;

}

void HydroMD::compute_drag(double dt,int iter) {
if (iter%couplingsteps==0) {
  for(int i=0;i<meso->size();i++) force[i]=vector3(0,0,0);//reset force 

  double i_sqrtdt=1.0/sqrt(dt*couplingsteps);
  double T = meso->T;

  receive();

  f = new vector3[ nohatoms ];//drag force
  double w[8];
  int cells[8];
  for(int i=0;i<nohatoms;i++) {
    interpolator(*meso, pos[i], w, cells);
    vector3 ivel(0,0,0);
    double W=0.0;
    for (int h=0;h<8;h++) {//interpolate hydro vel 
      ivel += w[h]*fc[ cells[h] ].P/fc[ cells[h] ].M;
      W += w[h];
    }
    double coeff = sqrt(2.0*BOLTZMAN*T*gamma)*i_sqrtdt;
    vector3 ff(normal(),normal(),normal());
    f[i] = - gamma*(vel[i] - ivel) + ff*coeff;//drag force on particle
    for (int h=0;h<8;h++)  force[ cells[h] ] += - w[h]/W*f[i]; //distribute momentum exchange over cells 
  }
  send();

  delete [] f;
}