forked from lijiext/lammps
add .cpp and .h to root src
This commit is contained in:
Ulf R. Pedersen
parent
75b3f34a58
commit
348c4eb7f3
|
|
@ -0,0 +1,245 @@
|
||||||
|
/*
|
||||||
|
fix_rhok.cpp
|
||||||
|
|
||||||
|
A fix to add harmonic potential that bias |rho(k)|.
|
||||||
|
|
||||||
|
The usage is as follows:
|
||||||
|
|
||||||
|
fix [name] [groupID] rhoK [nx] [ny] [nz] [kappa = spring constant] [rhoK0]
|
||||||
|
|
||||||
|
where k_i = (2 pi / L_i) * n_i
|
||||||
|
|
||||||
|
Written by Ulf Pedersen and Patrick Varilly, 4 Feb 2010
|
||||||
|
Tweaked for LAMMPS 15 Jan 2010 version by Ulf Pedersen, 19 Aug 2010
|
||||||
|
Tweaked again March 4th 2012 by Ulf R. Pedersen, and
|
||||||
|
September 2017 by Ulf R. Pedersen.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include "fix_rhok.h"
|
||||||
|
#include "error.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>
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
using namespace FixConst;
|
||||||
|
|
||||||
|
// Constructor: all the parameters to this fix specified in
|
||||||
|
// the LAMMPS input get passed in
|
||||||
|
FixRhok::FixRhok( LAMMPS* inLMP, int inArgc, char** inArgv )
|
||||||
|
: Fix( inLMP, inArgc, inArgv )
|
||||||
|
{
|
||||||
|
// Check arguments
|
||||||
|
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] );
|
||||||
|
}
|
||||||
|
|
||||||
|
FixRhok::~FixRhok()
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
// Methods that this fix implements
|
||||||
|
// --------------------------------
|
||||||
|
|
||||||
|
// Tells LAMMPS where this fix should act
|
||||||
|
int
|
||||||
|
FixRhok::setmask()
|
||||||
|
{
|
||||||
|
int mask = 0;
|
||||||
|
|
||||||
|
// This fix modifies forces...
|
||||||
|
mask |= POST_FORCE;
|
||||||
|
mask |= POST_FORCE_RESPA;
|
||||||
|
mask |= MIN_POST_FORCE;
|
||||||
|
|
||||||
|
// ...and potential energies
|
||||||
|
mask |= THERMO_ENERGY;
|
||||||
|
|
||||||
|
return mask;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*int FixRhok::setmask()
|
||||||
|
{
|
||||||
|
int mask = 0;
|
||||||
|
mask |= POST_FORCE;
|
||||||
|
mask |= POST_FORCE_RESPA;
|
||||||
|
mask |= MIN_POST_FORCE;
|
||||||
|
return mask;
|
||||||
|
}*/
|
||||||
|
|
||||||
|
|
||||||
|
// Initializes the fix at the beginning of a run
|
||||||
|
void
|
||||||
|
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 );
|
||||||
|
}
|
||||||
|
|
||||||
|
// Initial application of the fix to a system (when doing MD)
|
||||||
|
void
|
||||||
|
FixRhok::setup( int inVFlag )
|
||||||
|
{
|
||||||
|
if( strcmp( update->integrate_style, "verlet" ) == 0 )
|
||||||
|
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 );
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Initial application of the fix to a system (when doing minimization)
|
||||||
|
void
|
||||||
|
FixRhok::min_setup( int inVFlag )
|
||||||
|
{
|
||||||
|
post_force( inVFlag );
|
||||||
|
}
|
||||||
|
|
||||||
|
// Modify the forces calculated in the main force loop of ordinary MD
|
||||||
|
void
|
||||||
|
FixRhok::post_force( int inVFlag )
|
||||||
|
{
|
||||||
|
double **x = atom->x;
|
||||||
|
double **f = atom->f;
|
||||||
|
int *mask = atom->mask;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
|
// Loop over locally-owned atoms affected by this fix and calculate the
|
||||||
|
// partial rhoK's
|
||||||
|
mRhoKLocal[0] = 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
|
||||||
|
|
||||||
|
// 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];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Forces in RESPA loop
|
||||||
|
void
|
||||||
|
FixRhok::post_force_respa( int inVFlag, int inILevel, int inILoop )
|
||||||
|
{
|
||||||
|
if( inILevel == mNLevelsRESPA - 1 )
|
||||||
|
post_force( inVFlag );
|
||||||
|
}
|
||||||
|
|
||||||
|
// Forces in minimization loop
|
||||||
|
void
|
||||||
|
FixRhok::min_post_force( int inVFlag )
|
||||||
|
{
|
||||||
|
post_force( inVFlag );
|
||||||
|
}
|
||||||
|
|
||||||
|
// Compute the change in the potential energy induced by this fix
|
||||||
|
double
|
||||||
|
FixRhok::compute_scalar()
|
||||||
|
{
|
||||||
|
double rhoK = sqrt( mRhoKGlobal[0]*mRhoKGlobal[0]
|
||||||
|
+ mRhoKGlobal[1]*mRhoKGlobal[1] );
|
||||||
|
|
||||||
|
return 0.5 * mKappa * (rhoK - mRhoK0) * (rhoK - mRhoK0);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Compute the ith component of the vector
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,81 @@
|
||||||
|
/*
|
||||||
|
fix_rhok.h
|
||||||
|
|
||||||
|
A fix to do umbrella sampling on rho(k).
|
||||||
|
|
||||||
|
The usage is as follows:
|
||||||
|
|
||||||
|
fix [name] [groupID] rhoKUmbrella [kx] [ky] [kz] [kappa = spring constant] [rhoK0]
|
||||||
|
|
||||||
|
where k_i = (2 pi / L_i) * n_i
|
||||||
|
|
||||||
|
Written by Ulf Pedersen and Patrick Varilly, 4 Feb 2010
|
||||||
|
Tweaked for LAMMPS 15 Jan 2010 version by Ulf Pedersen, 19 Aug 2010
|
||||||
|
*/
|
||||||
|
|
||||||
|
#ifdef FIX_CLASS
|
||||||
|
|
||||||
|
FixStyle(rhok,FixRhok)
|
||||||
|
|
||||||
|
#else
|
||||||
|
|
||||||
|
#ifndef __FIX_RHOK__
|
||||||
|
#define __FIX_RHOK__
|
||||||
|
|
||||||
|
#include "fix.h"
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class FixRhok : public Fix
|
||||||
|
{
|
||||||
|
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();
|
||||||
|
|
||||||
|
// Methods that this fix implements
|
||||||
|
// --------------------------------
|
||||||
|
|
||||||
|
// Tells LAMMPS where this fix should act
|
||||||
|
int setmask();
|
||||||
|
|
||||||
|
// Initializes the fix at the beginning of a run
|
||||||
|
void init();
|
||||||
|
|
||||||
|
// Initial application of the fix to a system (when doing MD / minimization)
|
||||||
|
void setup( int inVFlag );
|
||||||
|
void min_setup( int inVFlag );
|
||||||
|
|
||||||
|
// Modify the forces calculated in the main force loop, either when
|
||||||
|
// doing usual MD, RESPA MD or minimization
|
||||||
|
void post_force( int inVFlag );
|
||||||
|
void post_force_respa( int inVFlag, int inILevel, int inILoop );
|
||||||
|
void min_post_force( int inVFlag );
|
||||||
|
|
||||||
|
// 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
|
||||||
|
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];
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace LAMMPS_NS
|
||||||
|
|
||||||
|
#endif // __FIX_RHOK__
|
||||||
|
#endif // FIX_CLASS
|
||||||
|
|
||||||
Loading…
Reference in New Issue