Skip to content

Commit

Permalink
Merge pull request #224 from jonaslandsgesell/master
Browse files Browse the repository at this point in the history 
Changed Metadynamics Algorithm, added relaxation steps
  • Loading branch information
@fweik
fweik committed Apr 17, 2015
2 parents 41d07a9 + f1cc4d3 commit e7b7971
Show file tree
Hide file tree
Showing 4 changed files with 34 additions and 19 deletions.
1 change: 1 addition & 0 deletions doc/ug/run.tex
View file
Expand Up @@ -455,6 +455,7 @@ \section{Metadynamics}
\item[\var{profile\_list}] Tcl list of a previous metadynamics
profile.
\item[\var{force\_list}] Tcl list of a previous metadynamics force.
\item[\var{numrelaxationsteps}] number of relaxation steps before setting a new hill.
\end{arguments}

\subsubsection{Details on usage}
Expand Down
15 changes: 10 additions & 5 deletions src/core/metadynamics.cpp
View file
Expand Up @@ -49,6 +49,8 @@ int meta_pid2 = -1;
double meta_bias_height = 0.001;
/** bias width */
double meta_bias_width = 0.5;
/** number of relaxation steps **/
int meta_num_relaxation_steps = -1;

/** REACTION COORDINATE */
/** RC min */
Expand Down Expand Up @@ -163,18 +165,21 @@ void meta_perform()
// reaction coordinate value
meta_val_xi = sqrt(sqrlen(meta_cur_xi));
// Update free energy profile and biased force
meta_acc_fprofile[i] -= calculate_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);
meta_acc_force[i] -= calculate_deriv_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);
if(int(sim_time/time_step)%meta_num_relaxation_steps==0){
meta_acc_fprofile[i] -= calculate_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);
meta_acc_force[i] -= calculate_deriv_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);
}

// direction of the bias force
unit_vector(meta_cur_xi,meta_apply_direction);
} else if (meta_switch == META_REL_Z) {
// reaction coordinate value: relative height of z_pid1 with respect to z_pid2
meta_val_xi = -1.*meta_cur_xi[2];
// Update free energy profile and biased force
meta_acc_fprofile[i] -= calculate_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);
meta_acc_force[i] -= calculate_deriv_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);

if(int(sim_time/time_step)%meta_num_relaxation_steps==0){
meta_acc_fprofile[i] -= calculate_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);
meta_acc_force[i] -= calculate_deriv_lucy(meta_xi_min+i*meta_xi_step,meta_val_xi);
}
// direction of the bias force (-1 to be consistent with META_DIST: from 1 to 2)
meta_apply_direction[0] = meta_apply_direction[1] = 0.;
meta_apply_direction[2] = -1.;
Expand Down
9 changes: 8 additions & 1 deletion src/core/metadynamics.hpp
View file
Expand Up @@ -53,9 +53,14 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#define META_REL_Z 2

/**********************************
* exported variables
* needed global variables
**********************************/
extern double sim_time;
extern double time_step;

/**********************************
* exported variables
**********************************/
/** Flag - turn metadynamics on */
extern int meta_switch;
/** pid of particle 1 */
Expand All @@ -66,6 +71,8 @@ extern int meta_pid2;
extern double meta_bias_height;
/** bias width */
extern double meta_bias_width;
/** number of relaxation steps before a new potential is set **/
extern int meta_num_relaxation_steps;

/** REACTION COORDINATE */
/** RC min */
Expand Down
28 changes: 15 additions & 13 deletions src/tcl/metadynamics_tcl.cpp
View file
Expand Up @@ -109,20 +109,20 @@ int tclcommand_metadynamics_parse_off(Tcl_Interp *interp, int argc, char **argv)

int tclcommand_metadynamics_parse_distance(Tcl_Interp *interp, int argc, char **argv)
{
int pid1, pid2, dbins;
int pid1, pid2, dbins, numrelaxationsteps;
double dmin, dmax, bheight, bwidth, fbound;

/* check number of arguments */
if (argc < 8) {
if (argc < 11) {
Tcl_AppendResult(interp, "wrong # args: should be \n\"",
argv[0]," ",argv[1]," <pid1> <pid2> <d_min> <d_max> <b_height> <b_width> <f_bound> <d_bins>\"", (char *)NULL);
argv[0]," ",argv[1]," <pid1> <pid2> <d_min> <d_max> <b_height> <b_width> <f_bound> <d_bins> <num_relaxation_steps>\"", (char *)NULL);
return (TCL_ERROR);
}

/* check argument types */
if ( !ARG_IS_I(2, pid1) || !ARG_IS_I(3, pid2) || !ARG_IS_D(4, dmin) || !ARG_IS_D(5, dmax) ||
!ARG_IS_D(6, bheight) || !ARG_IS_D(7, bwidth) || !ARG_IS_D(8, fbound) || !ARG_IS_I(9, dbins)) {
Tcl_AppendResult(interp, argv[0]," ",argv[1]," needs two INTS, five DOUBLES, and one INT", (char *)NULL);
!ARG_IS_D(6, bheight) || !ARG_IS_D(7, bwidth) || !ARG_IS_D(8, fbound) || !ARG_IS_I(9, dbins) || !ARG_IS_I(10, numrelaxationsteps ) ) {
Tcl_AppendResult(interp, argv[0]," ",argv[1]," needs two INTS, five DOUBLES, and two INTS in this order", (char *)NULL);
return (TCL_ERROR);
}

Expand All @@ -131,7 +131,7 @@ int tclcommand_metadynamics_parse_distance(Tcl_Interp *interp, int argc, char **
return (TCL_ERROR);
}

if (dmin < 0 || dmax < 0 || dmax < dmin || bheight < 0 || bwidth < 0 || fbound < 0 || dbins < 0) {
if (dmin < 0 || dmax < 0 || dmax < dmin || bheight < 0 || bwidth < 0 || fbound < 0 || dbins < 0 || numrelaxationsteps <0) {
Tcl_AppendResult(interp, "check parameters: inconcistency somewhere", (char *)NULL);
return (TCL_ERROR);
}
Expand All @@ -150,26 +150,27 @@ int tclcommand_metadynamics_parse_distance(Tcl_Interp *interp, int argc, char **
meta_xi_num_bins = dbins;

meta_switch = META_DIST;
meta_num_relaxation_steps = numrelaxationsteps;

return (TCL_OK);
}

int tclcommand_metadynamics_parse_relative_z(Tcl_Interp *interp, int argc, char **argv)
{
int pid1, pid2, dbins;
int pid1, pid2, dbins, numrelaxationsteps;
double dmin, dmax, bheight, bwidth, fbound;

/* check number of arguments */
if (argc < 8) {
if (argc < 11) {
Tcl_AppendResult(interp, "wrong # args: should be \n\"",
argv[0]," ",argv[1]," <pid1> <pid2> <z_min> <z_max> <b_height> <b_width> <f_bound> <z_bins>\"", (char *)NULL);
argv[0]," ",argv[1]," <pid1> <pid2> <z_min> <z_max> <b_height> <b_width> <f_bound> <z_bins> <num_relaxation_steps>\"", (char *)NULL);
return (TCL_ERROR);
}

/* check argument types */
if ( !ARG_IS_I(2, pid1) || !ARG_IS_I(3, pid2) || !ARG_IS_D(4, dmin) || !ARG_IS_D(5, dmax) ||
!ARG_IS_D(6, bheight) || !ARG_IS_D(7, bwidth) || !ARG_IS_D(8, fbound) || !ARG_IS_I(9, dbins)) {
Tcl_AppendResult(interp, argv[0]," ",argv[1]," needs two INTS, five DOUBLES, and one INT", (char *)NULL);
!ARG_IS_D(6, bheight) || !ARG_IS_D(7, bwidth) || !ARG_IS_D(8, fbound) || !ARG_IS_I(9, dbins) || !ARG_IS_I(10, numrelaxationsteps) ) {
Tcl_AppendResult(interp, argv[0]," ",argv[1]," needs two INTS, five DOUBLES, and two INTS in this order", (char *)NULL);
return (TCL_ERROR);
}

Expand All @@ -178,7 +179,7 @@ int tclcommand_metadynamics_parse_relative_z(Tcl_Interp *interp, int argc, char
return (TCL_ERROR);
}

if (dmax < dmin || bheight < 0 || bwidth < 0 || fbound < 0 || dbins < 0) {
if (dmax < dmin || bheight < 0 || bwidth < 0 || fbound < 0 || dbins < 0 || numrelaxationsteps <0) {
Tcl_AppendResult(interp, "check parameters: inconcistency somewhere", (char *)NULL);
return (TCL_ERROR);
}
Expand All @@ -197,6 +198,7 @@ int tclcommand_metadynamics_parse_relative_z(Tcl_Interp *interp, int argc, char
meta_xi_num_bins = dbins;

meta_switch = META_REL_Z;
meta_num_relaxation_steps = numrelaxationsteps;

return (TCL_OK);
}
Expand Down Expand Up @@ -268,7 +270,7 @@ int tclcommand_metadynamics_parse_load_stat(Tcl_Interp *interp, int argc, char *
Tcl_ResetResult(interp);
Tcl_SplitList(interp, argv[1], &tmp_argc, &tmp_argv);
realloc_doublelist(&profile, profile.n = tmp_argc);
printf("profile.n %d, meta_xi_num_bins %d\n",profile.n,meta_xi_num_bins);
//printf("profile.n %d, meta_xi_num_bins %d\n",profile.n,meta_xi_num_bins);
/* Now check that the number of items parsed is equal to the number of bins */
/* If there's one extra line, assume it's an empty line */
if (profile.n == meta_xi_num_bins+1)
Expand Down

0 comments on commit e7b7971

Please sign in to comment.