| 67 | 
  | 
    double Axy = hmat(0,0) * hmat(1, 1); | 
| 68 | 
  | 
    double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert); | 
| 69 | 
  | 
    double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert); | 
| 70 | 
< | 
    eta(0,0) -= Axy * (sx - surfaceTension) / (NkBT*tb2); | 
| 71 | 
< | 
    eta(1,1) -= Axy * (sy - surfaceTension) / (NkBT*tb2); | 
| 70 | 
> | 
    eta(0,0) -= dt2* Axy * (sx - surfaceTension) / (NkBT*tb2); | 
| 71 | 
> | 
    eta(1,1) -= dt2* Axy * (sy - surfaceTension) / (NkBT*tb2); | 
| 72 | 
  | 
    eta(2,2) += dt2 *  instaVol * (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); | 
| 73 | 
  | 
    oldEta = eta;   | 
| 74 | 
  | 
  } | 
| 80 | 
  | 
    prevEta = eta; | 
| 81 | 
  | 
    double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert); | 
| 82 | 
  | 
    double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert); | 
| 83 | 
< | 
    eta(0,0) = oldEta(0, 0) - Axy * (sx -surfaceTension) / (NkBT*tb2); | 
| 84 | 
< | 
    eta(1,1) = oldEta(1, 1) - Axy * (sy -surfaceTension) / (NkBT*tb2); | 
| 83 | 
> | 
    eta(0,0) = oldEta(0, 0) - dt2 * Axy * (sx -surfaceTension) / (NkBT*tb2); | 
| 84 | 
> | 
    eta(1,1) = oldEta(1, 1) - dt2 * Axy * (sy -surfaceTension) / (NkBT*tb2); | 
| 85 | 
  | 
    eta(2,2) = oldEta(2, 2) + dt2 *  instaVol * | 
| 86 | 
  | 
            (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); | 
| 87 | 
  | 
  } |