From 348c4eb7f30c441b1443aaeaf5317d74bd49e80b Mon Sep 17 00:00:00 2001
From: "Ulf R. Pedersen" <urp@nsm-i165>
Date: Thu, 28 Sep 2017 18:18:28 +0200
Subject: [PATCH] add .cpp and .h to root src
---
src/fix_rhok.cpp | 245 +++++++++++++++++++++++++++++++++++++++++++++++
src/fix_rhok.h | 81 ++++++++++++++++
2 files changed, 326 insertions(+)
create mode 100644 src/fix_rhok.cpp
create mode 100644 src/fix_rhok.h
diff --git a/src/fix_rhok.cpp b/src/fix_rhok.cpp
new file mode 100644
index 0000000000..6f5bd45a1d
--- /dev/null
+++ b/src/fix_rhok.cpp
@@ -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;
+}
diff --git a/src/fix_rhok.h b/src/fix_rhok.h
new file mode 100644
index 0000000000..3e0625430f
--- /dev/null
+++ b/src/fix_rhok.h
@@ -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
+
--
GitLab