jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make fix rhok more similar to other LAMMPS sources 

- re-indent to 2 blanks
- white space cleanup
- use force->numeric() and force->inumeric() instead of atof() and atoi()
- include system headers before local LAMMPS headers
This commit is contained in:
Axel Kohlmeyer 2017-10-03 10:10:38 -04:00
parent f07719e924
commit 251f28049a
2 changed files with 117 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

209
src/USER-MISC/fix_rhok.cpp
View File

@ -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;
}

30
src/USER-MISC/fix_rhok.h
View File

@ -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