From bc221879eab43829b10f553eb41e7ab6f014c132 Mon Sep 17 00:00:00 2001 From: Mikael Lund Date: Tue, 21 Feb 2017 23:16:50 +0100 Subject: [PATCH] Added vmd dipole moment update script for trajectories --- include/faunus/_xytable.h | 277 -------------------------------------- scripts/vmd-dipole.tcl | 44 ++++++ 2 files changed, 44 insertions(+), 277 deletions(-) delete mode 100644 include/faunus/_xytable.h create mode 100644 scripts/vmd-dipole.tcl diff --git a/include/faunus/_xytable.h b/include/faunus/_xytable.h deleted file mode 100644 index e0fb8bf23..000000000 --- a/include/faunus/_xytable.h +++ /dev/null @@ -1,277 +0,0 @@ -#ifndef FAU_XYTABLE_H -#define FAU_XYTABLE_H - -#ifndef SWIG -#include -#endif - -namespace Faunus { - /*! - * \brief Class for a 2D table, XY data etc. - * - * The types of X and Y are arbitrarily chosen.\n - * X values: - * - can be smaller than zero (see constructor) - * - are assumed equidistant - * - * The internal data vector will be automatically resized to fit - * the added data. However, for better memory utilization it is - * possible to specify a maximum x value. (if this is exceeded the - * vector will be resized anyway). - * - * Examples\n - * \code - * #include "average.h" - * #include "xytable.h" - * ... - * xytable f(0.1, -10); - * float x=-5.6; - * f(x)=10; - * cout << f(x); // --> 10 - * - * xytable > h(0.5); - * float r=7.5; - * h(r)+=2; - * h(r)+=4; - * cout << h(r).avg(); // --> 3 - * \endcode - * - * \author Mikael Lund - * \todo Check if the STL resize command assign zero to new data - */ - template - class xytable_center { - protected: - int x2i(TX x) { - int i=int( invxres*(x-xmin)+0.5); - if (i>=y.size()) - y.resize(i+1); - return i; - } - - inline int x2ifloor(TX x) { - return int( invxres*(x-xmin) ); - } - - public: - typedef typename std::vector::iterator yiter; - std::vector y; - TX xres, //!< Distance between x values - invxres, //!< Inverse distance between x values - xmin; //!< Minimum x value - - //! \param resolution Distance between x values - //! \param xminimum Minimum x value (can be smaller than zero) - //! \param xmaximum Maximum x value (for better memory optimization) - xytable_center(TX xminimum, TX xmaximum, TX resolution) { - init(resolution, xminimum, xmaximum); - } - - void init(TX xminimum, TX xmaximum, TX resolution) { - y.clear(); - xres=resolution; - invxres=1/xres; - xmin=xminimum; - if (xmaximum>0) - y.resize( int(std::abs(xmaximum-xminimum)/resolution) ); - } - - //! Returns the x value with the highest y value - TX x_atmax_y() { - yiter iter = std::max_element( y.begin(), y.end() ); - int i=iter-y.begin(); // convert iterator to index integer - return i*xres + xmin; - } - - //! Convenient data access - TY& operator()(TX val) { return y.at( x2i(val) ); } - - TX x(int i) { - return i*xres+xmin; - } - - /*! \brief Linear interpolation - * \warning If one of the reference values is Inf, nan is returned - */ - TY interpolate(TX val) { - int i = x2ifloor(val); - TY y_low = y[i]; - TY slope = (y[i+1]-y_low) * invxres; - return y_low + slope * (val-x(i)); - } - - //! Max x-value in set - TX xmax() { return (y.size()-1)*xres+xmin; } - - std::vector gety() { return y; } - - std::vector getx() { - std::vector v; - for (int i=0; i> n >> xmin >> xres; // extract number of points, xmin, xres - getline(f,s); // ignore aditional columns - y.resize(n); - for (int i=0; i> y[i]; - } - f.close(); - return true; - } - return false; - } - }; - - /*! - * \brief Class for a 2D table, XY data, etc. #2 - * Based on xytable_center with some modifications, - * here x-values correspond to the lower bound of the slice - * for which the value is stored instead of the middle - * of the slice - * - * \warning Does not automatically increase table size - * but points to the last slice instead - * \warning Function fails if x < xmin - * \author Chris Evers, Lund 2011-05-13 - */ - template - class xytable { - protected: - inline int x2i(TX x) { return floor( invxres*(x-xmin)); } - - public: - typedef typename std::vector::iterator yiter; - std::vector y; - TX xres, //!< Distance between x values - invxres, //!< Inverse distance between x values - xmin; //!< Minimum x value - int imax; //!< Index of maximum x value - - //! \param resolution Distance between x values - //! \param xminimum Minimum x value (can be smaller than zero) - //! \param xmaximum Maximum x value (for better memory optimization) - xytable(TX xminimum, TX xmaximum, TX resolution) { - init(resolution, xminimum, xmaximum); - } - - void init(TX xminimum, TX xmaximum, TX resolution) { - y.clear(); - xres=resolution; - invxres=1/resolution; - xmin=xminimum; - imax=ceil(std::abs(xmaximum-xminimum)/resolution)-1; - if (xmaximum>0) - y.resize(imax+1); // last one is used if x=xmaximum - } - - //! Returns the x value with the highest y value - TX x_atmax_y() { - yiter iter = std::max_element( y.begin(), y.end() ); - int i=iter-y.begin(); // convert iterator to index integer - return i*xres + xmin; - } - - //! Convenient data access - TY& operator()(TX val) { return y.at( x2i(val) ); } - - TX x(int i) { - return i*xres+xmin; - } - - /*! \brief Linear interpolation - * \warning If one of the reference values is Inf, nan is returned - */ - TY interpolate(TX val) { - int i = x2i(val); - TY y_low = y[i]; - TY slope = (y[i+1]-y_low) * invxres; - return y_low + slope * (val-x(i)); - } - - //! Max x-value in set - TX xmax() { return imax*xres+xmin; } - - std::vector gety() { return y; } - - std::vector getx() { - std::vector v; - for (int i=0; i> n >> xmin >> xres; - y.resize(n); - for (auto y_i : y) - f >> y_i; - f.close(); - return true; - } - return false; - } - }; -}//namespace -#endif diff --git a/scripts/vmd-dipole.tcl b/scripts/vmd-dipole.tcl new file mode 100644 index 000000000..cf64a2eb3 --- /dev/null +++ b/scripts/vmd-dipole.tcl @@ -0,0 +1,44 @@ +# Script to draw dipole moment vectors for resid's +# and update for every frame in the tracjectory. +# Dipole moment based on the geometric center. +# Mikael Lund, September 2009. +# +# Usage: +# Specify resid range in for-loop below. +# $ source dipol.tcl +# $ enabletrace + +proc vmd_draw_arrow {mol start end} { + set end [vecadd $start [vecscale 0.4 $end]] + set middle [vecadd $start [vecscale 0.8 [vecsub $end $start]]] + set start [vecadd $start [vecscale -0.9 [vecsub $end $start]]] + graphics $mol cylinder $start $middle radius 1.4 + #graphics $mol cylinder $start $end radius 1.4 + graphics $mol cone $middle $end radius 2.2 +} + +proc dipmom {args} { + graphics 0 delete all + for { set i 1 } { $i <= 1 } { incr i } { + set sel1 [atomselect top "resid $i"] + set muvec [measure dipole $sel1] + set cm [measure center $sel1] + #set vec [veclength $muvec] + graphics 0 color red + vmd_draw_arrow top $cm $muvec + } +} + +proc enabletrace {} { + global vmd_frame + #trace add variable ::vmd_frame([molinfo top]) write dipmom + trace variable vmd_frame([molinfo top]) w dipmom +} + +proc disabletrace {} { + global vmd_frame + #trace remove variable ::vmd_frame([molinfo top]) write dipmom + trace vdelete vmd_frame([molinfo top]) w dipmom +} + +