--- trunk/OOPSE/libmdtools/ExtendedSystem.cpp 2003/04/04 01:57:11 454 +++ trunk/OOPSE/libmdtools/ExtendedSystem.cpp 2003/04/04 19:16:11 457 @@ -14,7 +14,8 @@ ExtendedSystem::ExtendedSystem( SimInfo &info ) { atoms = info.atoms; nMols = info.n_mol; molecules = info.molecules; - zeta = 0; + zeta = 0.0; + epsilonDot = 0.0; } @@ -22,36 +23,36 @@ void ExtendedSystem::NoseHooverNVT( double dt ){ } -void ExtendedSystem::NoseHooverNVT( double dt ){ +void ExtendedSystem::NoseHooverNVT( double dt, double ke ){ // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 int i; - double kB, keconverter, NkBT, zetaScale, ke_temp; + double NkBT, zetaScale, ke_temp; double vx, vy, vz, jx, jy, jz; - - kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K - keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2 / K - - ke_temp = getKinetic() * keconverter; + const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K + const double e_convert = 4.184e-4; // to convert ke from kcal/mol to + // amu*Ang^2*fs^-2/K + + ke_temp = ke * e_convert; NkBT = (double)getNDF() * kB * targetTemp; - // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin & + // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin // qmass is set in the parameter file - zeta += dt*((ke_temp*2 - NkBT)/qmass); + + zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); zetaScale = zeta * dt; // perform thermostat scaling on linear velocities and angular momentum - for(i = 0; i < n_atoms; i++){ vx = atoms[i]->get_vx(); vy = atoms[i]->get_vy(); vz = atoms[i]->get_vz(); - atoms[i]->set_vx(vx - zetaScale * vx); - atoms[i]->set_vy(vy - zetaScale * vy); - atoms[i]->set_vz(vz - zetaScale * vz); + atoms[i]->set_vx(vx * (1.0 - zetaScale)); + atoms[i]->set_vy(vy * (1.0 - zetaScale)); + atoms[i]->set_vz(vz * (1.0 - zetaScale)); } if( n_oriented ){ @@ -65,65 +66,63 @@ void ExtendedSystem::NoseHooverNVT( double dt ){ jy = dAtom->getJy(); jz = dAtom->getJz(); - dAtom->setJx( jx - zetaScale * jx); - dAtom->setJy( jy - zetaScale * jy); - dAtom->setJz( jz - zetaScale * jz); + dAtom->setJx(jx * (1.0 - zetaScale)); + dAtom->setJy(jy * (1.0 - zetaScale)); + dAtom->setJz(jz * (1.0 - zetaScale)); } } } } -void ExtendedSystem::NoseHooverAndersonNPT(double pressure, double ke, - double dt, double temp ) { +void ExtendedSystem::NoseHooverAndersonNPT( double dt, + double ke, + double p_int ) { // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500 - int i, j, degrees_freedom; - double pressure, dt, temp, pressure_units, epsilonScale; - double ke, kB, vxi, vyi, vzi, pressure_ext; - double boxx_old, boxy_old, boxz_old; - double keconverter, NkBT, zetaScale, ke_temp; - double jxi, jyi, jzi, scale; + double oldBox[3]; + double newBox[3]; + const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K + const double p_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1 + const double e_convert = 4.184e-4; // to convert ke from kcal/mol to + // amu*Ang^2*fs^-2/K - kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K - pressure_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1 - degrees_freedom = 6*nmol; // number of degrees of freedom for the system - keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2/K + double p_ext; - pressure_ext = pressure * pressure_units; - volume = boxx*boxy*boxz; - ke_temp = ke * keconverter; - NkBT = degrees_freedom*kB*temp; + p_ext = targetPressure * p_units; + p_mol = p_int * p_units; + getBox(oldBox); + + volume = oldBox[0]*oldBox[1]*oldBox[2]; + + ke_temp = ke * e_convert; + NkBT = (double)getNDF() * kB * targetTemp; + // propogate the strain rate - epsilon_dot += dt * ( (p_mol - pressure_ext)*volume - / (tau_relax*tau_relax * kB * targetTemp) ); + epsilonDot += dt * ((p_mol - p_ext) * volume / + (tauRelax*tauRelax * kB * targetTemp) ); // determine the change in cell volume - scale = pow( (1.0 + dt * 3.0 * epsilon_dot), (1.0 / 3.0)); + scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0)); - volume = volume * pow(scale, 3.0); + newBox[0] = oldBox[0] * scale; + newBox[1] = oldBox[1] * scale; + newBox[2] = oldBox[2] * scale; + volume = newBox[0]*newBox[1]*newBox[2]; // perform affine transform to update positions with volume fluctuations - affine_transform( scale ); + this->AffineTransform( oldBox, newBox ); - // save old lengths and update box size - boxx_old = boxx; - boxy_old = boxy; - boxz_old = boxz; - - boxx = boxx_old*scale; - boxy = boxy_old*scale; - boxz = boxz_old*scale; - epsilonScale = epsilonDot * dt; // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin // qmass is set in the parameter file - zeta += dt * ( (ke_temp*2 - NkBT) / qmass ); + + zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); zetaScale = zeta * dt; // apply barostating and thermostating to velocities and angular momenta @@ -133,9 +132,9 @@ void ExtendedSystem::NoseHooverAndersonNPT(double pres vy = atoms[i]->get_vy(); vz = atoms[i]->get_vz(); - atoms[i]->set_vx(vx * (1.0 - zetaScale * epsilonScale)); - atoms[i]->set_vy(vy * (1.0 - zetaScale * epsilonScale)); - atoms[i]->set_vz(vz * (1.0 - zetaScale * epsilonScale)); + atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale)); + atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale)); + atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale)); } if( n_oriented ){ @@ -157,41 +156,47 @@ void ExtendedSystem::AffineTransform( double scale ){ } } -void ExtendedSystem::AffineTransform( double scale ){ +void ExtendedSystem::AffineTransform( double oldBox[3], double newBox[3] ){ int i; - double boxx_old, boxy_old, boxz_old, percentScale; - double boxx_num, boxy_num, boxz_num, rxi, ryi, rzi; - double[3] r; + double r[3]; + double boxNum[3]; + double percentScale[3]; + double rxi, ryi, rzi; // first determine the scaling factor from the box size change - percentScale = (boxx - boxx_old)/boxx_old; + percentScale[0] = (newBox[0] - oldBox[0]) / oldBox[0]; + percentScale[1] = (newBox[1] - oldBox[1]) / oldBox[1]; + percentScale[2] = (newBox[2] - oldBox[2]) / oldBox[2]; - for (i=0; i < nMols; i++) { molecules[i]->getCOM(r); - // find the minimum image coordinates of the molecular centers of mass: + // find the minimum image coordinates of the molecular centers of mass: - - boxx_num = boxx_old*copysign(1.0,r[0])*(double)(int)(fabs(r[0]/boxx_old)+0.5); + boxNum[0] = oldBox[0] * copysign(1.0,r[0]) * + (double)(int)(fabs(r[0]/oldBox[0]) + 0.5); - boxx_num = boxx_old*dsign(1.0d0,rx(i))*int(abs(rx(i)/boxx_old)+0.5d0); - boxy_num = boxy_old*dsign(1.0d0,ry(i))*int(abs(ry(i)/boxy_old)+0.5d0); - boxz_num = boxz_old*dsign(1.0d0,rz(i))*int(abs(rz(i)/boxz_old)+0.5d0); + boxNum[1] = oldBox[1] * copysign(1.0,r[1]) * + (double)(int)(fabs(r[1]/oldBox[1]) + 0.5); - rxi = rx(i) - boxx_num; - ryi = ry(i) - boxy_num; - rzi = rz(i) - boxz_num; + boxNum[2] = oldBox[2] * copysign(1.0,r[2]) * + (double)(int)(fabs(r[2]/oldBox[2]) + 0.5); + rxi = r[0] - boxNum[0]; + ryi = r[1] - boxNum[1]; + rzi = r[2] - boxNum[2]; + // update the minimum image coordinates using the scaling factor - rxi = rxi + rxi*percentScale; - ryi = ryi + ryi*percentScale; - rzi = rzi + rzi*percentScale; + rxi += rxi*percentScale[0]; + ryi += ryi*percentScale[1]; + rzi += rzi*percentScale[2]; - rx(i) = rxi + boxx_num; - ry(i) = ryi + boxy_num; - rz(i) = rzi + boxz_num; + r[0] = rxi + boxNum[0]; + r[1] = ryi + boxNum[1]; + r[2] = rzi + boxNum[2]; + + molecules[i]->moveCOM(r); } }