diff --git a/doc/ug/run.tex b/doc/ug/run.tex index 52ffaff336a..193a0d6d364 100644 --- a/doc/ug/run.tex +++ b/doc/ug/run.tex @@ -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} diff --git a/src/core/metadynamics.cpp b/src/core/metadynamics.cpp index cdf8d20e2af..06d055e36e4 100644 --- a/src/core/metadynamics.cpp +++ b/src/core/metadynamics.cpp @@ -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 */ @@ -163,8 +165,10 @@ 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); @@ -172,9 +176,10 @@ void meta_perform() // 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.; diff --git a/src/core/metadynamics.hpp b/src/core/metadynamics.hpp index 83cbadc3aa9..48b2a3aab20 100644 --- a/src/core/metadynamics.hpp +++ b/src/core/metadynamics.hpp @@ -53,9 +53,14 @@ along with this program. If not, see . #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 */ @@ -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 */ diff --git a/src/tcl/metadynamics_tcl.cpp b/src/tcl/metadynamics_tcl.cpp index 8122928245e..0a3562eec6f 100644 --- a/src/tcl/metadynamics_tcl.cpp +++ b/src/tcl/metadynamics_tcl.cpp @@ -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]," \"", (char *)NULL); + argv[0]," ",argv[1]," \"", (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); } @@ -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); } @@ -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]," \"", (char *)NULL); + argv[0]," ",argv[1]," \"", (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); } @@ -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); } @@ -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); } @@ -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)