From 7cebec382db9401b15b8046bc8dffadb51f1f9db Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Sun, 13 Sep 2009 08:13:13 -0400 Subject: [PATCH] Added piston tension model. Also typographic cleanup 'assert (' -> 'assert('. Currently broken: stiffness calculation with the tension model. --- src/sawsim.nw | 169 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 157 insertions(+), 12 deletions(-) diff --git a/src/sawsim.nw b/src/sawsim.nw index 9917130..c247d7c 100644 --- a/src/sawsim.nw +++ b/src/sawsim.nw @@ -264,6 +264,7 @@ following models have been implemented: \item Hookean spring (Appendix \ref{sec.hooke}), \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}). + \item Piston (Appendix \ref{sec.piston_tension}), \end{packed_item} The tension model and parameters used for a particular domain depend @@ -418,7 +419,7 @@ you accidentally assign incompatible domains to the same group you should get an appropriate error message. <>= -#define NUM_TENSION_MODELS 5 +#define NUM_TENSION_MODELS 6 @ @@ -806,7 +807,8 @@ tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = { <>, <>, <>, - <> + <>, + <> }; @ @@ -947,7 +949,7 @@ void get_options(int argc, char **argv, extern char *optarg; extern int optind, optopt, opterr; - assert (argc > 0); + assert(argc > 0); #ifdef DEBUG for (i=0; i 0); + assert(stiffness >= 0); return stiffness; } @ @@ -3695,6 +3697,7 @@ The interface is defined in [[tension_model.h]], the implementation in [[tension <> <> <> +<> <> <> #endif /* TENSION_MODEL_H */ @@ -3764,7 +3767,7 @@ double const_tension_handler(double x, void *pdata) list_t *list = data->group_tension_params; double F; - assert (x >= 0.0); + assert(x >= 0.0); list = head(list); assert(list != NULL); /* empty active group?! */ F = ((const_tension_param_t *)list->d)->F; @@ -3829,7 +3832,7 @@ char const_tension_param_string[] = "0"; @ <>= -{&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t} +{&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", num_const_tension_params, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t} @ \subsection{Hooke} @@ -3896,7 +3899,7 @@ double hooke_handler(double x, void *pdata) { hooke_param_t param = {0}; - assert (x >= 0.0); + assert(x >= 0.0); hooke_param_grouper((tension_handler_data_t *)pdata, ¶m); return param.k*x; } @@ -3909,7 +3912,7 @@ double inverse_hooke_handler(double F, void *pdata) { hooke_param_t param = {0}; - assert (F >= 0.0); + assert(F >= 0.0); hooke_param_grouper((tension_handler_data_t *)pdata, ¶m); return F/param.k; } @@ -4132,6 +4135,7 @@ double wlc_handler(double x, void *pdata) tension_handler_data_t *data = (tension_handler_data_t *)pdata; wlc_param_t param = {0}; + assert(x >= 0.0); wlc_param_grouper(data, ¶m); return wlc(x, data->env->T, param.p, param.L); } @@ -4214,7 +4218,7 @@ char wlc_param_string[]="0.39e-9,28.4e-9"; @ <>= -{wlc_handler, inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t} +{&wlc_handler, &inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t} @ Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696. @@ -4348,6 +4352,145 @@ char fjc_param_string[]="0.5e-9,200"; {fjc_handler, NULL, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t} @ +\subsection{Piston} +\label{sec.piston_tension} + +An infinitely stretchable domain with no tension for extensions less +than a particular threshold $L$, and infinite tension for $x>L$. The +tension at $x=L$ is undefined, but may be any positive value. The +model is called the ``piston'' model because it resembles a +frictionless piston sliding in a rigid cylinder of length $L$. + +Because the tension at $x=L$ is undefined, the user must make sure +domains with this tension model are never the initial active group, +because it would break [[tension_balance()]] in [[find_tension()]] +(see Section \ref{sec.tension_balance}). + +<>= +<> +<> +<> +<> +@ + +<>= +<> +<> +<> +<> + +@ + +First we define a function that computes the effective piston parameters +(stored in a single [[piston_tension_param_t]]) for the entire group. +<>= +static void piston_tension_param_grouper(tension_handler_data_t *data, + piston_tension_param_t *param) { + list_t *list = data->group_tension_params; + double L=0; + + list = head(list); + assert(list != NULL); /* empty active group?! */ + while (list != NULL) { + L += ((piston_tension_param_t *)list->d)->L; + list = list->next; + } + param->L = L; +} + +<>= +extern double piston_tension_handler(double x, void *pdata); +@ +<>= +double piston_tension_handler(double x, void *pdata) +{ + tension_handler_data_t *data = (tension_handler_data_t *)pdata; + piston_tension_param_t param = {0}; + double F; + + piston_tension_param_grouper(data, ¶m); + if (x <= param.L) F = 0; + else F = HUGE_VAL; +#ifdef DEBUG + fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F); +#endif + return F; +} + +@ + +<>= +extern double inverse_piston_tension_handler(double x, void *pdata); +@ +<>= +double inverse_piston_tension_handler(double F, void *pdata) +{ + tension_handler_data_t *data = (tension_handler_data_t *)pdata; + piston_tension_param_t param = {0}; + + assert(F >= 0.0); + piston_tension_param_grouper(data, ¶m); + if (F == 0.0) return 0; + return param.L; +} + +@ + +<>= +typedef struct piston_tension_param_struct { + double L; /* length of piston in m */ +} piston_tension_param_t; + +@ + +<>= +extern void *string_create_piston_tension_param_t(char **param_strings); +extern void destroy_piston_tension_param_t(void *p); + +@ + +<>= +piston_tension_param_t *create_piston_tension_param_t(double L) +{ + piston_tension_param_t *ret + = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t)); + assert(ret != NULL); + ret->L = L; + return ret; +} + +void *string_create_piston_tension_param_t(char **param_strings) +{ + return create_piston_tension_param_t( + safe_strtod(param_strings[0],__FUNCTION__)); +} + +void destroy_piston_tension_param_t(void *p) +{ + if (p) + free(p); +} + +@ + +<>= +extern int num_piston_tension_params; +extern char *piston_tension_param_descriptions[]; +extern char piston_tension_param_string[]; + +@ + +<>= +int num_piston_tension_params = 1; +char *piston_tension_param_descriptions[] = {"length L, m"}; +char piston_tension_param_string[] = "0"; + +@ + +<>= +{&piston_tension_handler, &inverse_piston_tension_handler, "piston", "a domain that slides without tension for xL.", num_piston_tension_params, piston_tension_param_descriptions, piston_tension_param_string, &string_create_piston_tension_param_t, &destroy_piston_tension_param_t} +@ + \subsection{Tension model implementation} <>= @@ -4376,6 +4519,7 @@ char fjc_param_string[]="0.5e-9,200"; <> <> <> +<> @ <>= @@ -4384,6 +4528,7 @@ char fjc_param_string[]="0.5e-9,200"; <> <> <> +<> @ <>= @@ -4392,9 +4537,9 @@ char fjc_param_string[]="0.5e-9,200"; <> <> <> +<> @ - \subsection{Utilities} The tension models can get complicated, and users may want to reassure @@ -4545,7 +4690,7 @@ void get_options(int argc, char **argv, environment_t *env, extern char *optarg; extern int optind, optopt, opterr; - assert (argc > 0); + assert(argc > 0); /* setup defaults */ @@ -5894,7 +6039,7 @@ void get_options(int argc, char **argv, environment_t *env, extern char *optarg; extern int optind, optopt, opterr; - assert (argc > 0); + assert(argc > 0); /* setup defaults */ -- 2.26.2