diff --git a/src/USER-PINNING/README b/src/USER-PINNING/README new file mode 100644 index 0000000000..9f158b1884 --- /dev/null +++ b/src/USER-PINNING/README @@ -0,0 +1,13 @@ +This bias potential is used to study solid-liquid transitions with the interface pinning method. + +Reference: + +[Ulf R. Pedersen, J. Chem. Phys. 139, 104102 (2013)] + +Usage snip: + + fix [fix-name] [groupID] rhoKUmbrella [nx] [ny] [nz] [kappa] [anchor-point] + thermo_style custom step temp pzz pe lz f_umbrella f_umbrella[1] f_umbrella[2] f_umbrella[3] + +where the parameters set the harmonic bias potential U=0.5*kappa*(|rho_k|-anchor-point)^2 +with the wave-vector elements of rho_k to k_x = (2 pi / L_x) * n_x, k_y = (2 pi / L_y) * n_y and k_z = (2 pi / L_z) * n_z. diff --git a/src/USER-PINNING/fix_rhoKUmbrella.cpp b/src/USER-PINNING/fix_rhoKUmbrella.cpp new file mode 100644 index 0000000000..e62eecb31c --- /dev/null +++ b/src/USER-PINNING/fix_rhoKUmbrella.cpp @@ -0,0 +1,246 @@ +/* + fix_rhoK_umbrella.cpp + + A fix to do umbrella sampling on rho(k). + + The usage is as follows: + + fix [name] [groupID] rhoKUmbrella [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 Pedersen. + */ + +#include "fix_rhoKUmbrella.h" +#include "error.h" +#include "update.h" +#include "respa.h" +#include "atom.h" +#include "domain.h" + +#include +#include +#include +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +// Constructor: all the parameters to this fix specified in +// the LAMMPS input get passed in +FixRhoKUmbrella::FixRhoKUmbrella( 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] ); +} + +FixRhoKUmbrella::~FixRhoKUmbrella() +{ +} + +// Methods that this fix implements +// -------------------------------- + +// Tells LAMMPS where this fix should act +int +FixRhoKUmbrella::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 FixRhoKUmbrella::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 +FixRhoKUmbrella::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 +FixRhoKUmbrella::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 +FixRhoKUmbrella::min_setup( int inVFlag ) +{ + post_force( inVFlag ); +} + +// Modify the forces calculated in the main force loop of ordinary MD +void +FixRhoKUmbrella::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 ); + + // WARNING!!!!! < \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 + // is intensive. + // + // Don't forget this two years from now when you change the system size!!! + 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 +FixRhoKUmbrella::post_force_respa( int inVFlag, int inILevel, int inILoop ) +{ + if( inILevel == mNLevelsRESPA - 1 ) + post_force( inVFlag ); +} + +// Forces in minimization loop +void +FixRhoKUmbrella::min_post_force( int inVFlag ) +{ + post_force( inVFlag ); +} + +// Compute the change in the potential energy induced by this fix +double +FixRhoKUmbrella::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 +FixRhoKUmbrella::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; +} diff --git a/src/USER-PINNING/fix_rhoKUmbrella.h b/src/USER-PINNING/fix_rhoKUmbrella.h new file mode 100644 index 0000000000..d549d31d94 --- /dev/null +++ b/src/USER-PINNING/fix_rhoKUmbrella.h @@ -0,0 +1,81 @@ +/* + fix_rhoK_umbrella.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(rhoKUmbrella,FixRhoKUmbrella) + +#else + +#ifndef __FIX_RHOKUMBRELLA__ +#define __FIX_RHOKUMBRELLA__ + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixRhoKUmbrella : public Fix +{ +public: + // Constructor: all the parameters to this fix specified in + // the LAMMPS input get passed in + FixRhoKUmbrella( LAMMPS* inLMP, int inArgc, char** inArgv ); + virtual ~FixRhoKUmbrella(); + + // 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_RHOKUMBRELLA__ +#endif // FIX_CLASS +