Commit 3495141d authored by Tomáš Trnka's avatar Tomáš Trnka
Browse files

Fix the target kinetic energy of the NH barostat 

The cell momenta should be thermostatted to kT per barostat degree
of freedom (d^2 in general, d*(d-1) without rotations), according to
Shinoda et al. 2004 (doi:10.1103/PhysRevB.69.134103) Eqn. 1 and Martyna,
Tobias, Klein (JCP 1994, doi:10.1063/1.467468 section II.D).
parent d0ba8e1d

src/fix_nh.cpp

+23 −9
Original line number Diff line number Diff line
@@ -1446,7 +1446,7 @@ double FixNH::compute_scalar() 
  double volume;

  double energy;

  double kt = boltz * t_target;

  double lkt_press = kt;

  double lkt_press = 0.0;

  int ich;

  if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;

  else volume = domain->xprd * domain->yprd;
@@ -1477,15 +1477,21 @@ double FixNH::compute_scalar() 
  //       sum is over barostatted dimensions



  if (pstat_flag) {

    for (i = 0; i < 3; i++)

      if (p_flag[i])

    for (i = 0; i < 3; i++) {

      if (p_flag[i]) {

        energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] +

          p_hydro*(volume-vol0) / (pdim*nktv2p);

        lkt_press += kt;

      }

    }



    if (pstyle == TRICLINIC) {

      for (i = 3; i < 6; i++)

        if (p_flag[i])

      for (i = 3; i < 6; i++) {

        if (p_flag[i]) {

          energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i];

          lkt_press += kt;

        }

      }

    }



    // extra contributions from thermostat chain for barostat
@@ -1818,10 +1824,10 @@ void FixNH::nhc_temp_integrate() 


void FixNH::nhc_press_integrate()

{

  int ich,i;

  int ich,i,pdof;

  double expfac,factor_etap,kecurrent;

  double kt = boltz * t_target;

  double lkt_press = kt;

  double lkt_press;



  // Update masses, to preserve initial freq, if flag set
@@ -1850,14 +1856,22 @@ void FixNH::nhc_press_integrate() 
  }



  kecurrent = 0.0;

  pdof = 0;

  for (i = 0; i < 3; i++)

    if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];

    if (p_flag[i]) {

      kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];

      pdof++;

    }



  if (pstyle == TRICLINIC) {

    for (i = 3; i < 6; i++)

      if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];

      if (p_flag[i]) {

        kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];

        pdof++;

      }

  }



  lkt_press = pdof * kt;

  etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];



  double ncfac = 1.0/nc_pchain;

CRA Git | Maintained and supported by SUSTech CRA and CCSE