From 60c67b07dc3785f3db79a417c1680796589abc0f Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Wed, 26 Jul 2017 10:45:11 -0400
Subject: [PATCH] import updated fix msst file with some additional cleanup and
simplification
---
src/SHOCK/fix_msst.cpp | 91 +++++++++++++++++++++++-------------------
1 file changed, 51 insertions(+), 40 deletions(-)
diff --git a/src/SHOCK/fix_msst.cpp b/src/SHOCK/fix_msst.cpp
index b66f0ed4d3..66a648cd17 100644
--- a/src/SHOCK/fix_msst.cpp
+++ b/src/SHOCK/fix_msst.cpp
@@ -446,7 +446,7 @@ void FixMSST::initial_integrate(int vflag)
{
int i,k;
double p_msst; // MSST driving pressure
- double vol,TS,TS_term,escale_term;
+ double vol;
int nlocal = atom->nlocal;
int *mask = atom->mask;
@@ -469,12 +469,16 @@ void FixMSST::initial_integrate(int vflag)
// must convert energy to mv^2 units
if (dftb) {
- double TS_dftb = fix_external->compute_vector(0);
- TS = force->ftm2v*TS_dftb;
+ const double TS_dftb = fix_external->compute_vector(0);
+ const double TS = force->ftm2v*TS_dftb;
+ // update S_elec terms and compute TS_dot via finite differences
+ S_elec_2 = S_elec_1;
+ S_elec_1 = S_elec;
+ const double Temp = temperature->compute_scalar();
+ S_elec = TS/Temp;
+ TS_dot = Temp*(3.0*S_elec-4.0*S_elec_1+S_elec_2)/(2.0*update->dt);
+ TS_int += (update->dt*TS_dot);
if (update->ntimestep == 1) T0S0 = TS;
- } else {
- TS = 0.0;
- T0S0 = 0.0;
}
// compute new pressure and volume
@@ -484,16 +488,6 @@ void FixMSST::initial_integrate(int vflag)
couple();
vol = compute_vol();
- // update S_elec terms and compute TS_dot via finite differences
-
- S_elec_2 = S_elec_1;
- S_elec_1 = S_elec;
- double Temp = temperature->compute_scalar();
- S_elec = TS/Temp;
- TS_dot = Temp*(3.0*S_elec-4.0*S_elec_1+S_elec_2)/(2.0*update->dt);
- TS_int += (update->dt*TS_dot);
- //TS_int += (update->dt*TS_dot)/total_mass;
-
// compute etot + extra terms for conserved quantity
double e_scale = compute_etotal() + compute_scalar();
@@ -530,9 +524,9 @@ void FixMSST::initial_integrate(int vflag)
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
- double C = f[i][k] * force->ftm2v / mass[type[i]];
- TS_term = TS_dot/(mass[type[i]]*velocity_sum);
- escale_term = force->ftm2v*beta*(e0-e_scale) /
+ const double C = f[i][k] * force->ftm2v / mass[type[i]];
+ const double TS_term = TS_dot/(mass[type[i]]*velocity_sum);
+ const double escale_term = force->ftm2v*beta*(e0-e_scale) /
(mass[type[i]]*velocity_sum);
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
@@ -540,7 +534,7 @@ void FixMSST::initial_integrate(int vflag)
old_velocity[i][k] = v[i][k];
if ( k == direction ) D -= 2.0 * omega[sd] / vol;
if ( fabs(dthalf * D) > 1.0e-06 ) {
- double expd = exp(D * dthalf);
+ const double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
@@ -553,15 +547,15 @@ void FixMSST::initial_integrate(int vflag)
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
- double C = f[i][k] * force->ftm2v / mass[type[i]];
+ const double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
old_velocity[i][k] = v[i][k];
if ( k == direction ) {
- D = D - 2.0 * omega[sd] / vol;
+ D -= 2.0 * omega[sd] / vol;
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
- double expd = exp(D * dthalf);
+ const double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
@@ -590,16 +584,16 @@ void FixMSST::initial_integrate(int vflag)
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
- double C = f[i][k] * force->ftm2v / mass[type[i]];
- TS_term = TS_dot/(mass[type[i]]*velocity_sum);
- escale_term = force->ftm2v*beta*(e0-e_scale) /
+ const double C = f[i][k] * force->ftm2v / mass[type[i]];
+ const double TS_term = TS_dot/(mass[type[i]]*velocity_sum);
+ const double escale_term = force->ftm2v*beta*(e0-e_scale) /
(mass[type[i]]*velocity_sum);
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
D += escale_term - TS_term;
if ( k == direction ) D -= 2.0 * omega[sd] / vol;
if ( fabs(dthalf * D) > 1.0e-06 ) {
- double expd = exp(D * dthalf);
+ const double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
@@ -612,14 +606,14 @@ void FixMSST::initial_integrate(int vflag)
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
- double C = f[i][k] * force->ftm2v / mass[type[i]];
+ const double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
if ( k == direction ) {
- D = D - 2.0 * omega[sd] / vol;
+ D -= 2.0 * omega[sd] / vol;
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
- double expd = exp(D * dthalf);
+ const double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
@@ -669,7 +663,6 @@ void FixMSST::final_integrate()
{
int i;
double p_msst; // MSST driving pressure
- double TS_term,escale_term;
// v update only for atoms in MSST group
@@ -687,22 +680,38 @@ void FixMSST::final_integrate()
double e_scale = compute_etotal() + compute_scalar();
+ // for DFTB, extract TS_dftb from fix external
+ // must convert energy to mv^2 units
+
+ if (dftb) {
+ const double TS_dftb = fix_external->compute_vector(0);
+ const double TS = force->ftm2v*TS_dftb;
+ S_elec_2 = S_elec_1;
+ S_elec_1 = S_elec;
+ const double Temp = temperature->compute_scalar();
+ // update S_elec terms and compute TS_dot via finite differences
+ S_elec = TS/Temp;
+ TS_dot = Temp*(3.0*S_elec-4.0*S_elec_1+S_elec_2)/(2.0*update->dt);
+ TS_int += (update->dt*TS_dot);
+ if (update->ntimestep == 1) T0S0 = TS;
+ }
+
// propagate particle velocities 1/2 step
if (dftb) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( int k = 0; k < 3; k++ ) {
- double C = f[i][k] * force->ftm2v / mass[type[i]];
- TS_term = TS_dot/(mass[type[i]]*velocity_sum);
- escale_term = force->ftm2v*beta*(e0-e_scale) /
+ const double C = f[i][k] * force->ftm2v / mass[type[i]];
+ const double TS_term = TS_dot/(mass[type[i]]*velocity_sum);
+ const double escale_term = force->ftm2v*beta*(e0-e_scale) /
(mass[type[i]]*velocity_sum);
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
D += escale_term - TS_term;
if ( k == direction ) D -= 2.0 * omega[sd] / vol;
if ( fabs(dthalf * D) > 1.0e-06 ) {
- double expd = exp(D * dthalf);
+ const double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
@@ -715,14 +724,14 @@ void FixMSST::final_integrate()
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( int k = 0; k < 3; k++ ) {
- double C = f[i][k] * force->ftm2v / mass[type[i]];
+ const double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
if ( k == direction ) {
- D = D - 2.0 * omega[sd] / vol;
+ D -= 2.0 * omega[sd] / vol;
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
- double expd = exp(D * dthalf);
+ const double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
@@ -748,7 +757,7 @@ void FixMSST::final_integrate()
( v0 - vol )/( v0 * v0 );
double A = total_mass * ( p_current[sd] - p0 - p_msst ) /
( qmass * nktv2p * mvv2e );
- double B = total_mass * mu / ( qmass * vol );
+ const double B = total_mass * mu / ( qmass * vol );
// prevent blow-up of the volume
@@ -950,7 +959,9 @@ double FixMSST::compute_scalar()
// subtract off precomputed TS_int integral value
- energy -= TS_int;
+ if (dftb) { // TS_int == 0 for non DFTB calculations
+ energy -= TS_int;
+ }
return energy;
}
--
GitLab