diff --git a/src/USER-MISC/fix_rhok.cpp b/src/USER-MISC/fix_rhok.cpp
index df8b8bd6678818b08a8c6e62ee8c3e4b9219045d..58b0e95a97e0abe10ac7d48e1e93807a0c5a16ab 100644
--- a/src/USER-MISC/fix_rhok.cpp
+++ b/src/USER-MISC/fix_rhok.cpp
@@ -13,24 +13,24 @@
Contributing author: Ulf R. Pedersen, ulf@urp.dk
------------------------------------------------------------------------- */
-#include "fix_rhok.h"
-#include "error.h"
-#include "citeme.h"
-#include "update.h"
-#include "respa.h"
-#include "atom.h"
-#include "domain.h"
-
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include <assert.h>
#include <math.h>
+#include "fix_rhok.h"
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "respa.h"
+#include "update.h"
+#include "citeme.h"
+
using namespace LAMMPS_NS;
using namespace FixConst;
-static const char cite_flow_gauss[] =
+static const char cite_fix_rhok[] =
"Bias on the collective density field (fix rhok):\n\n"
"@Article{pedersen_jcp139_104102_2013,\n"
"title = {Direct calculation of the solid-liquid Gibbs free energy difference in a single equilibrium simulation},\n"
@@ -48,40 +48,35 @@ static const char cite_flow_gauss[] =
FixRhok::FixRhok( LAMMPS* inLMP, int inArgc, char** inArgv )
: Fix( inLMP, inArgc, inArgv )
{
-
- if (lmp->citeme) lmp->citeme->add(cite_flow_gauss);
-
+
+ if (lmp->citeme) lmp->citeme->add(cite_fix_rhok);
+
// Check arguments
- if( inArgc != 8 )
- error->all(FLERR,"Illegal fix rhoKUmbrella command" );
+ if( inArgc != 8 )
+ error->all(FLERR,"Illegal fix rhoKUmbrella command" );
// Set up fix flags
scalar_flag = 1; // have compute_scalar
vector_flag = 1; // have compute_vector...
size_vector = 3; // ...with this many components
global_freq = 1; // whose value can be computed at every timestep
- //scalar_vector_freq = 1; // OLD lammps: whose value can be computed at every timestep
thermo_energy = 1; // this fix changes system's potential energy
extscalar = 0; // but the deltaPE might not scale with # of atoms
extvector = 0; // neither do the components of the vector
// Parse fix options
- int n[3];
-
- n[0] = atoi( inArgv[3] );
- n[1] = atoi( inArgv[4] );
- n[2] = atoi( inArgv[5] );
-
- mK[0] = n[0]*(2*M_PI / (domain->boxhi[0] - domain->boxlo[0]));
- mK[1] = n[1]*(2*M_PI / (domain->boxhi[1] - domain->boxlo[1]));
- mK[2] = n[2]*(2*M_PI / (domain->boxhi[2] - domain->boxlo[2]));
-
- mKappa = atof( inArgv[6] );
- mRhoK0 = atof( inArgv[7] );
-}
+ int n[3];
-FixRhok::~FixRhok()
-{
+ n[0] = force->inumeric(FLERR,inArgv[3]);
+ n[1] = force->inumeric(FLERR,inArgv[4]);
+ n[2] = force->inumeric(FLERR,inArgv[5]);
+
+ mK[0] = n[0]*(2*M_PI / (domain->boxhi[0] - domain->boxlo[0]));
+ mK[1] = n[1]*(2*M_PI / (domain->boxhi[1] - domain->boxlo[1]));
+ mK[2] = n[2]*(2*M_PI / (domain->boxhi[2] - domain->boxlo[2]));
+
+ mKappa = force->numeric(FLERR,inArgv[6]);
+ mRhoK0 = force->numeric(FLERR,inArgv[7]);
}
// Methods that this fix implements
@@ -110,20 +105,20 @@ FixRhok::init()
{
// RESPA boilerplate
if( strcmp( update->integrate_style, "respa" ) == 0 )
- mNLevelsRESPA = ((Respa *) update->integrate)->nlevels;
-
- // Count the number of affected particles
- int nThisLocal = 0;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
- for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
- if( mask[i] & groupbit ) { // ...only those affected by this fix
- nThisLocal++;
- }
- }
- MPI_Allreduce( &nThisLocal, &mNThis,
- 1, MPI_INT, MPI_SUM, world );
- mSqrtNThis = sqrt( mNThis );
+ mNLevelsRESPA = ((Respa *) update->integrate)->nlevels;
+
+ // Count the number of affected particles
+ int nThisLocal = 0;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
+ if( mask[i] & groupbit ) { // ...only those affected by this fix
+ nThisLocal++;
+ }
+ }
+ MPI_Allreduce( &nThisLocal, &mNThis,
+ 1, MPI_INT, MPI_SUM, world );
+ mSqrtNThis = sqrt( mNThis );
}
// Initial application of the fix to a system (when doing MD)
@@ -131,13 +126,13 @@ void
FixRhok::setup( int inVFlag )
{
if( strcmp( update->integrate_style, "verlet" ) == 0 )
- post_force( inVFlag );
+ post_force( inVFlag );
else
- {
- ((Respa *) update->integrate)->copy_flevel_f( mNLevelsRESPA - 1 );
- post_force_respa( inVFlag, mNLevelsRESPA - 1,0 );
- ((Respa *) update->integrate)->copy_f_flevel( mNLevelsRESPA - 1 );
- }
+ {
+ ((Respa *) update->integrate)->copy_flevel_f( mNLevelsRESPA - 1 );
+ post_force_respa( inVFlag, mNLevelsRESPA - 1,0 );
+ ((Respa *) update->integrate)->copy_f_flevel( mNLevelsRESPA - 1 );
+ }
}
// Initial application of the fix to a system (when doing minimization)
@@ -159,54 +154,54 @@ FixRhok::post_force( int inVFlag )
// Loop over locally-owned atoms affected by this fix and calculate the
// partial rhoK's
mRhoKLocal[0] = 0.0;
- mRhoKLocal[1] = 0.0;
+ mRhoKLocal[1] = 0.0;
for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
- if( mask[i] & groupbit ) { // ...only those affected by this fix
+ if( mask[i] & groupbit ) { // ...only those affected by this fix
- // rho_k = sum_i exp( - i k.r_i )
- mRhoKLocal[0] += cos( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
- mRhoKLocal[1] -= sin( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
- }
+ // rho_k = sum_i exp( - i k.r_i )
+ mRhoKLocal[0] += cos( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
+ mRhoKLocal[1] -= sin( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
+ }
+ }
+
+ // Now calculate mRhoKGlobal
+ MPI_Allreduce( mRhoKLocal, mRhoKGlobal,
+ 2, MPI_DOUBLE, MPI_SUM, world );
+
+ // Info: < \sum_{i,j} e^{-ik.(r_i - r_j)} > ~ N, so
+ // we define rho_k as (1 / sqrt(N)) \sum_i e^{-i k.r_i}, so that
+ // <rho_k^2> is intensive.
+ mRhoKGlobal[0] /= mSqrtNThis;
+ mRhoKGlobal[1] /= mSqrtNThis;
+
+ // We'll need magnitude of rho_k
+ double rhoK = sqrt( mRhoKGlobal[0]*mRhoKGlobal[0]
+ + mRhoKGlobal[1]*mRhoKGlobal[1] );
+
+ for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
+ if( mask[i] & groupbit ) { // ...only those affected by this fix
+
+ // Calculate forces
+ // U = kappa/2 ( |rho_k| - rho_k^0 )^2
+ // f_i = -grad_i U = -kappa ( |rho_k| - rho_k^0 ) grad_i |rho_k|
+ // grad_i |rho_k| = Re( rho_k* (-i k e^{-i k . r_i} / sqrt(N)) ) / |rho_k|
+ //
+ // In terms of real and imag parts of rho_k,
+ //
+ // Re( rho_k* (-i k e^{-i k . r_i}) ) =
+ // (- Re[rho_k] * sin( k . r_i ) - Im[rho_k] * cos( k . r_i )) * k
+
+ double sinKRi = sin( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
+ double cosKRi = cos( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
+
+ double prefactor = mKappa * ( rhoK - mRhoK0 ) / rhoK
+ * (-mRhoKGlobal[0]*sinKRi - mRhoKGlobal[1]*cosKRi) / mSqrtNThis;
+ f[i][0] -= prefactor * mK[0];
+ f[i][1] -= prefactor * mK[1];
+ f[i][2] -= prefactor * mK[2];
+ }
}
-
- // Now calculate mRhoKGlobal
- MPI_Allreduce( mRhoKLocal, mRhoKGlobal,
- 2, MPI_DOUBLE, MPI_SUM, world );
-
- // Info: < \sum_{i,j} e^{-ik.(r_i - r_j)} > ~ N, so
- // we define rho_k as (1 / sqrt(N)) \sum_i e^{-i k.r_i}, so that
- // <rho_k^2> is intensive.
- mRhoKGlobal[0] /= mSqrtNThis;
- mRhoKGlobal[1] /= mSqrtNThis;
-
- // We'll need magnitude of rho_k
- double rhoK = sqrt( mRhoKGlobal[0]*mRhoKGlobal[0]
- + mRhoKGlobal[1]*mRhoKGlobal[1] );
-
- for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
- if( mask[i] & groupbit ) { // ...only those affected by this fix
-
- // Calculate forces
- // U = kappa/2 ( |rho_k| - rho_k^0 )^2
- // f_i = -grad_i U = -kappa ( |rho_k| - rho_k^0 ) grad_i |rho_k|
- // grad_i |rho_k| = Re( rho_k* (-i k e^{-i k . r_i} / sqrt(N)) ) / |rho_k|
- //
- // In terms of real and imag parts of rho_k,
- //
- // Re( rho_k* (-i k e^{-i k . r_i}) ) =
- // (- Re[rho_k] * sin( k . r_i ) - Im[rho_k] * cos( k . r_i )) * k
-
- double sinKRi = sin( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
- double cosKRi = cos( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
-
- double prefactor = mKappa * ( rhoK - mRhoK0 ) / rhoK
- * (-mRhoKGlobal[0]*sinKRi - mRhoKGlobal[1]*cosKRi) / mSqrtNThis;
- f[i][0] -= prefactor * mK[0];
- f[i][1] -= prefactor * mK[1];
- f[i][2] -= prefactor * mK[2];
- }
- }
}
// Forces in RESPA loop
@@ -214,7 +209,7 @@ void
FixRhok::post_force_respa( int inVFlag, int inILevel, int inILoop )
{
if( inILevel == mNLevelsRESPA - 1 )
- post_force( inVFlag );
+ post_force( inVFlag );
}
// Forces in minimization loop
@@ -228,9 +223,9 @@ FixRhok::min_post_force( int inVFlag )
double
FixRhok::compute_scalar()
{
- double rhoK = sqrt( mRhoKGlobal[0]*mRhoKGlobal[0]
- + mRhoKGlobal[1]*mRhoKGlobal[1] );
-
+ double rhoK = sqrt( mRhoKGlobal[0]*mRhoKGlobal[0]
+ + mRhoKGlobal[1]*mRhoKGlobal[1] );
+
return 0.5 * mKappa * (rhoK - mRhoK0) * (rhoK - mRhoK0);
}
@@ -238,13 +233,13 @@ FixRhok::compute_scalar()
double
FixRhok::compute_vector( int inI )
{
- if( inI == 0 )
- return mRhoKGlobal[0]; // Real part
- else if( inI == 1 )
- return mRhoKGlobal[1]; // Imagniary part
- else if( inI == 2 )
- return sqrt( mRhoKGlobal[0]*mRhoKGlobal[0]
- + mRhoKGlobal[1]*mRhoKGlobal[1] );
- else
- return 12345.0;
+ if( inI == 0 )
+ return mRhoKGlobal[0]; // Real part
+ else if( inI == 1 )
+ return mRhoKGlobal[1]; // Imagniary part
+ else if( inI == 2 )
+ return sqrt( mRhoKGlobal[0]*mRhoKGlobal[0]
+ + mRhoKGlobal[1]*mRhoKGlobal[1] );
+ else
+ return 12345.0;
}
diff --git a/src/USER-MISC/fix_rhok.h b/src/USER-MISC/fix_rhok.h
index 145b3e7c1e2cfc32119c0fa3ac83b538cfdb22c6..c950c08b1d14aa3b22f14ce98e9b5c65eb513750 100644
--- a/src/USER-MISC/fix_rhok.h
+++ b/src/USER-MISC/fix_rhok.h
@@ -15,8 +15,8 @@ FixStyle(rhok,FixRhok)
#else
-#ifndef __FIX_RHOK__
-#define __FIX_RHOK__
+#ifndef LMP_FIX_RHOK_H
+#define LMP_FIX_RHOK_H
#include "fix.h"
@@ -28,8 +28,8 @@ public:
// Constructor: all the parameters to this fix specified in
// the LAMMPS input get passed in
FixRhok( LAMMPS* inLMP, int inArgc, char** inArgv );
- virtual ~FixRhok();
-
+ virtual ~FixRhok() {};
+
// Methods that this fix implements
// --------------------------------
@@ -51,23 +51,23 @@ public:
// Compute the change in the potential energy induced by this fix
double compute_scalar();
-
- // Compute the ith component of the vector associated with this fix
+
+ // Compute the ith component of the vector associated with this fix
double compute_vector( int inI );
-
+
private:
// RESPA boilerplate
int mNLevelsRESPA;
// Defining parameters for this umbrella
- double mK[3], mKappa, mRhoK0;
-
- // Number of particles affected by the fix
- int mNThis;
- double mSqrtNThis;
-
- // Real and imaginary parts of rho_k := sum_i exp( - i k . r_i )
- double mRhoKLocal[2], mRhoKGlobal[2];
+ double mK[3], mKappa, mRhoK0;
+
+ // Number of particles affected by the fix
+ int mNThis;
+ double mSqrtNThis;
+
+ // Real and imaginary parts of rho_k := sum_i exp( - i k . r_i )
+ double mRhoKLocal[2], mRhoKGlobal[2];
};
} // namespace LAMMPS_NS