Added USER-PINNING package

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.

/*

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

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

	// <rho_k^2> 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;

}

/*

  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