\documentclass[a4paper]{spie}
\usepackage[]{graphicx}
%%%\usepackage{float}

%%===========================
%%========== Title ==========
%%===========================
\title{Finite element simulations of the bending of the IPMC sheet}

%%=============================
%%========== Authors ==========
%%=============================
\author{D. Pugal, H. Kasem\"agi, M. Kruusmaa, A. Aabloo, K. J. Kim
\skiplinehalf
Institute of Technology, Tartu University, Nooruse 1, 50411 Tartu, Estonia\\
University of Nevada, Reno
}
\authorinfo{Further author information: (Send correspondence to Alvo Aabloo)\\
Alvo Aabloo: Email: alvo.aabloo@ut.ee}

\pagestyle{plain}
\begin{document}
\maketitle

%%==============================
%%========== Abstract ==========
%%==============================

\begin{abstract}
Finite Element method is used to build and simulate an IPMC sheet. The physical bending of the realistic Nafion sheet due to the drift of counter-ions (e.g Na+) and dragged water in applied electric field are simulated. The other effects like  viscous resistance and diffusion of water, also the porous nature of the environment inside the Nafion are taken into account. The effect of the concentration of water and counter-ions near the electrodes are tried to relate to the physical bending of the IPMC sheet. Some additional effects like voltage drop at more distance points from contacts instead of constant electric field and electrolysis for more complex model are discussed.
\end{abstract}

\keywords{Electroactive polymers, EAP, Finite element method, Electrochemical-mechanical analysis, Actuator, Coupled problem}

%%==================================
%%========== Introduction ==========
%%==================================

\section{Introduction}

EAP-based electromechanical actuators are valuable for use in a number of applications starting with miniature robotics upto military and space. These actuators have light weigth, noiseless motion, simple mechanical constrution; large controlled displacement and good damage tolerance along with an ability to perform different movements like bending and contractions makes possible to use them as artificial muscles.

%Polymers whose shape can be purposely changed by applying some kind of electric signal to them are known as Electroactive Polymers (EAPs). They can be used as actuators or sensors. Large controlled displacements and good damage tolerance are achievable, similar to those of biological muscles, therefore such materials are also known as artificial muscles and they have possible application in the field of robotics \cite{cohen}.

%EAP-based electromechanical actuators have light weight, noiseless motion and simple mechanical construction while able to perform different movements like contraction, bending etc., and can be valuable for use in miniature robotics, medical, space, military and other applications where inflicting similar movement by the means of traditional electromechanic components would be too awkward.

%%%EAP generally utilise one of two mechanisms:
%%%\begin{itemize}
%%%	\item The movements of electronic EAP are caused by electrostatic forces. They are easy to keep in a steady position, but as a downside they require electric fields over 150 V/m to activate.
%%%	\item The shape of ionic EAP are changed by ion movements inside them, activated by considerably lower voltage of 1--5 V, but require continuous current to stay in a chosen position \cite{cohen}.
%%%\end{itemize}
%%%Our research group focuses on ionic polymer EAP materials, specifically ionic polymer metal composites (IPMC). Experiments were done with an underwater robot having two pectoral fins made of 0.2$\ldots$0.5 mm thick, 50 mm long and 3 to 12 mm wide electroactive polymer sheet covered with platinum coating (IPMC) provided by Musclesheet\texttrademark \cite{punn01}. The test confirmed the ability of the fins to generate thrust and move the robot forward. The results were found to depend largely on the properties of the IPMC.
%%%
%%%Each of the two 110-mm fins of our next prototype robot consisted 8 bottle-shaped muscles from IPMC sheet by Musclesheet\texttrademark: 40 mm long, 4 to 13 mm wide and 0.2$\ldots$0.5 mm thick \cite{punn02}. The fins were able to generate undulating motion and propel the body forward, although the speed was lower than that of ray-like fishes due to the smaller muscle mass and total mass ratio.
%%%
%%%We introduced a novel open-loop control method for IPMC artificial muscles to reduce their power consumption. The new method was found to be consistent with theoretical findings, tested experimentally and verified to the existing methods. As a result, the energy consumption of the muscles can be reduced upto 3.7 times and peak current consumption up to 2.6 times \cite{punn03}.

%In order to improve the properties of the IPMC artificial muscles, the surface resistance of IPMC actuators and sensors was studied \cite{punn04}. The results show surface resistance high correlation to bending curvature. The equivalent circuit model was suggested to describe the actual behaviour of IPMC sensors and actuators and the effects of the surface resistance were discussed.

%The next step from here is to incorporate molecular level computer simulations to better understand the behaviour of IPMCs. We have a lot of experiences of simulations of different polymeric system \cite{THOMAS-SSI-2001-143-83-87,THOMAS-JMATCHEM-2003-13-214-218,THOMAS-SSI-2004-168-249-254,THOMAS-CHEMMAT-2005-17-3673-3680,THOMAS-JMATCHEM-2005-15-1422-1428,THOMAS-JMATCHEM-2005-15-4338-4345}, mostly for lithium-ion polymer battery applications. Also we have parametrised force fields \cite{THOMAS-JCOMPHYS-2005-26-716-724} with Quantum Chemistry (QC) methods \cite{jensen} between different ions and ionic complexes for Molecular Dynamics (MD) simulations \cite{allen,rapaport,ercolessi,ohno,taopang}.

%In the present work a new program for generating initial configurations, i. e. atom/ion positions in the MD simulation box with Monte Carlo method is introduced. The properties of the program are closely described, the optimal values for program parameters are found and the performance of the program is analysed.
 
%%========================================
%%========== Simulation details ==========
%%========================================

\section{SIMULATION DETAILS}

%Monte Carlo Method is a method often used to study the macroscopic properties, e. g. phase transformation temperatures of structures consisting of many chaotically moving particles and numerical integration of intricite mathematical functions. Large number of random data points are created and a statistical analysis is used to find the macroscopic parameters. Data points can be uniformly distributed and taken into account with different weighs (simple sampling technique) or already generated with the weighs used as generation probability for increased efficiency (importance sampling technique) if the desired distribution function is quite far from uniform distribution. \cite{ohno} \textit{mcgen} falls to the latter category for although the variables of adding another monomer unit are generated with uniform distribution, a good unit configuration is used for long time unchanged while adding more units to it, but bad configurations will soon be replaced.

%The suitable probability function $p_{x}$ for simulation of a system of atoms is the Boltzmann distribution $p_{x}\sim\exp{\frac{-E_x}{kT}}$ where $x$ is the selected configuration, $E_{x}$ its potential energy, $T$ is the simulation temperature and $k$ is the Boltzmann constant \cite{leach}.

%When a long-chained macromolecul is randomly placed into a simulation cell with limited measures, the molecule probably overlaps itself in the beginning configuration, a situation that should be avoided \cite{leach}.

%Incremental adding of atoms to an open end of the chain in the means of monomers is similar to a radical chain polymerization process \cite{stuart} and would be better suited for macromolecules than random placement of the whole chain, so \textit{mcgen} adopts this kind of algorithm.

%This modified method comes together with change of energy $E_{k}$ caused by the change of the positions of atoms and their configuration potential energies, and with additional change of energy $E_{1}$due to the change of the number of atoms in the simulation cell. After each growth step a measurement will be taken of the total energy change $\Delta E=E_{1}+E_{k}$. 

%It has been proven \cite{METROPOLIS-JCHEMPHYS-21-1953-1087-1092} that the Metropolis method that does not change the number of atoms and therefore ensures that $E_{1}=0$ gives natural distribution of configurations with different potentials.

%When in our method $E_{1}<0$ as will often happen in case of creation of many covalent bonds with negative potential into an empty simulation box, the configurations distribution becomes less specific as all configurations with energies below the threshold $E_{k}\le E_{1}$ are always accepted.

%When the potential energy is biased so that $E_{1}>0$, then the configurations distribution function becomes abnormally narrow as the possibility of having $\Delta E<0$ is decreased. Also the probability function is multiplied by factor $0<\exp\frac{-E_{1}}{kT}<1$ thus decreasing the general unit acception rate and unnecessarily increasing simulation time.

%Therefore the model should be made in such a way that the value of $E_{1}$ be as close to zero as possible. The gradual growing methodology does not allow this value to be made exactly zero by simple means because the average energy increase per added monomer unit does depend on the current amount of free space in the cell by increasing the part of long-range energy, but some level of correction by adding an empirically derived constant is useful. The uneven distribution of dihedral angles along the polymer chain can be eliminated by following relaxation in Molecular Dynamics simulation \cite{colbourn}.

%%=====================================
%%========== program details ==========
%%=====================================

%\section{PROGRAM DETAILS}

%\subsection{General Properties}

%The current status of the program include the following features:

%\begin{itemize}
%	\item growth of single or multiple chains with the same and/or different structure;
%	\item growth of branched chains with constant or random distance between the branches along the backbone;
%	\item insertion of atoms and/or ions to the simulation box prior to the chain growth;
%	\item exclusion of preinserted atoms from the simulation box when they stay on the way of the growing chain;
%	\item geometrical constraints on the way of the growing chain; the chain can grow inside or outside of a sphere; stay on the one side of a wall (or walls) perpendicular to the coordinate axis (axes); stay inside or outside a rectangular box (each coordinate axis has a wall perpendicular to it); multiple constraints of the same or different type can be specified at the same time.
%\end{itemize}

%\subsection{Control Parameters}

%The growth of the chain is controlled by the following parameters:

%\begin{description}
%	\item[MAXTRYUNIT] aka. MTU - the number of unaccepted configurations for the current unit before the configuration of the previous accepted unit becomes unaccepted; the growth continues with trying the first unaccepted unit again;
%	\item[RTRY] - the number of unaccepted configurations for the current unit before the configurations of more than one last accepted units become unaccepted; the growth continues with trying the first unaccepted unit again;
%	\item[RSTEP] - the number of units to be unaccepted in the case of fixed takeback of units or the highest number of units taken back in the case of random takeback;
%	\item[RFIX] - takeback strategy; it determines if the number of takeback units is random or fixed.
%\end{description}

%\begin{figure}
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.75\textwidth]{../figures/mcgen-algorithm-060119.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:algorithm} Program algorithm}
%\end{figure}

%\subsection{The Algorithm}

%The algorithm is presented on Fig. . At first, the program reads structure and preinserted atom coordinates from input files. The growth of each molecule starts by creating the data structure in the memory. If the molecule is not immediately finished, the coordinates of the next unit are generated. At first, the unit is not accepted if the coordinates of any of the atoms in the unit fall into the region restricted by geometric constraints or coordinates overlap to any other atom within the sum of those two atoms' radiuses already present in the simulation box. Here, exceptional are the atoms connected by any type of bonds. When this test is passed, the energy difference $\Delta E$ between the new configuration and the previous one is calculated. If $\Delta E \le 0$, the unit in its new configuration is accepted at once. If $\Delta E > 0$, a random number $q \in (0,1)$ is compared to the probability $p = exp\frac{-\Delta E}{kT}$. The configuration is accepted in the case $q < p$ \cite{leach}.

%If one unit gets not accepted up to MAXTRYUNIT times sequentially in the algorithm, it is removed and a new configuration is generated for the last accepted unit. Each unit has another counter RTRY with non-zero initial value and every fulfilled MAXTRYUNIT is subtracted from the RTRY value of according unit. If this variable is counted down to zero for a unit, it is reset with its initial value and a number of units are taken back. This number is determined by RSTEP and RFIX values. RFIX defines the takeback strategy that can be either fixed or random. In the case of fixed takeback RSTEP is the constant number of units to take back. In the case of random strategy, the number of takeback units is chosen randomly from 2 upto the value of RSTEP.

%%============================================
%%========== Results and discussion ==========
%%============================================

\section{RESULTS AND DISCUSSION}

%\subsection{Model Parameters and Test Setup}

%Poly(ethylene oxide) (PEO) \cite{TADOKORO-MACMOL-6-1973-672-675} is chosen for test cases. PEO is a widely used prototype electrolyte host material for rechargeable lithium-ion polymer batteries and therefore well-studied both experimentally and theoretically. Test cases try to generate the amorphous phase of the polymer, since it is well-known that this phase has better ionic conductivity over the crystalline one; and this is one of the properties important for battery applications \cite{bib:gray,WRIGHT-ELCHIMACTA-43-1998-1137-1143,ARMAND-SOLSTATION-63-1994-309-319}.

%At present work, three types of simulation boxes are studied:
%\begin{itemize}
%	\item cubic simulation box with side length of 35 {\AA} and periodical boundary conditions at all three direction; this simulation box should contain a chain of 582 monomers;
%	\item cubic simulation box with side length of 35 {\AA} and periodical boundary conditions in two directions; the third dimension has following geometrical constraint: $0 \le x \le 17.5$ {\AA}, so the chain is free to grow between 2 walls in y- and z-direction; the target for this simulation box is to create a chain of 290 monomers as needed for a slab;
%	\item cubic simulation box with side length of 35 {\AA} and periodical boundary conditions in one direction; the second and third dimensions have following geometrical constraints: $0 \le x \le 17.5$ {\AA} and $0 \le y \le 17.5$ {\AA}, so the chain is free to grow between 4 walls in z-direction; the target for this simulation box is to create a chain of 290 monomers as needed for a wire.
%\end{itemize}
%The methyl groups are scheduled to add for chain ends. The force field for PEO was developed earlier \cite{THOMAS-JCHEMPHYS-101-1994-10064-10073}.

%The purpose of the test cases is to find optimal values of program control parameters MTU, RTRY, RSTEP and RFIX. Optimal values of control parameters should ensure the least computing expenses on generating the simulation box of desired density.

%The program receives 72 hours of real time to finish the chain or enough time to perform at least 50 million tries to add new monomer for each set of control parameters. In the first test series the value of MTU is fixed at 500, RTRY is varied in the range of 2000\ldots 16000 (1000\ldots 16000 for 3D periodicity and random takeback) and RSTEP in the range of 2\ldots 10. Then the best result of the first series are studied by varying MTU in the range of 100\ldots 900. In the analysis of the results the second series are compared and efficiencies calculated against the results for MTU=500.

%\subsection{Three-Dimensional Freedom}

%\subsubsection{Fixed takeback at MTU=500}

%The first testing step with MTU=500 and fixed takeback algorithm gives following results 
%(Fig.~\ref{fig:3Dfixtop})
%. In the region $RTRY<1000\cdot(RSTEP+1)$ takeback is probably too large-scale since the production of units decreases on increasing the proportion of takeback (decreased RTRYS and/or increased RSTEP). There is a possibility to abandon potentially good moves in this case.

%The effect on the result of RTRY changes does decrease on increasing RTRY value. Apparently the good moves are kept and less time spent to resolve hopeless moves. Those result on increasing RSTEP shows that takeback is overacting. Thus there is no point in using larger RSTEP values than 4 in the case of fixed takeback.

%%%========== fig:3Dfixtop ==========
%\begin{figure}[!ht]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/peak-unit-p3Dp0105293175r01-f-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:3Dfixtop}
%Relative number of generated monomers with fixed takeback at different RSTEP values and MTU=500 for 3-dimensional freedom}
%\end{figure}
%%%========== fig:3Drandomtop ==========
%\begin{figure}[!hb]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/peak-unit-p3Dp0105293175r01-r-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:3Drandomtop}
%Relative number of generated monomers with random takeback at different RSTEP values and MTU=500 for 3-dimensional freedom}
%\end{figure}

%%%========== fig:3DvariedMTU ==========
%\begin{figure}[h]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/analyys_yle_mtu_3D-RF-unit-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:3DvariedMTU}
%Relative number of generated monomers at selected RTRY/RSTEP values for 3-dimensional freedom}
%\end{figure}

%RSTEP=2 is the most effective for several RTRYS values, but these are accompanied by the most ineffective results. So short takeback does not ensure exiting from the 4-monomer traps.

%The best control parameters are RTRY=12000/RSTEP=3 (554 monomers or 95 \% of the maximum number of monomers), RTRY=13000/RSTEP=4 (543 monomers or 93 \%) and\\
%RTRY=7000/RSTEP=2 (562 monomers or 96 \%)
%%\\[5mm]
%%According to the minimal energy the maximum shifts to RTRY=14000 at RSTEP=4. The best RSTEP=5 is at RTRY=13000. This is somewhat surprising, indicating that program can resolve 5-monomer trap but on the expense of higher potential energy.

%\subsubsection{Random takeback at MTU=500}

%In the case of random takeback (Fig.) the region of hyperactive takeback ends at lower RTRY values (usually 3000\ldots 4000) than for fixed takeback. Appears that the region of premature takeback reaches to the end approximately at the same average number of abandoned units both for random and fixed cases. RSTEP=2 and RSTEP=3 show unstable behaviour. The best results in the broad range of RTRY gives RSTEP=5 and RSTEP=6 with average number of abandoned units of 3.5 and 4 respectively.

%The global maximum according to the longest chain is at RTRY=16000/RSTEP=4 (559 mo\-no\-mers) but this is located at the edge of the examined RTRY range and neighbouring RTRY values have worse results. The next best is at RTRY=8000/RSTEP=6 (549 monomers or 94 \%) and this is extremely stable against the variations of RSTEP; RSTEP$\pm$2 and RTRY$\pm$2000 do not decrease the result below 90 \%.

%One can conclude that random takeback gives about 1 \% better results than fixed takeback.

%\subsubsection{Varied MTU}
%The selected parameters (Fig.) are RTRY=2000/RSTEP=2, RTRY=7000/RSTEP=2, RTRY=12000/\\RSTEP=3 and RTRY=13000/RSTEP=4 for fixed takeback; RTRY=10000/RSTEP=5, RTRY=12000/\\RSTEP=5, RTRY=13000/RSTEP=6 and RTRY=15000/RSTEP=8 for random takeback.

%In the case of fixed takeback, the dependence of the results on the varied MTU is weak. The change is less than 15 \% in the whole RTRY range. For RTRY=2000/RSTEP=2 and RTRY=7000/RSTEP=2 the varying of MTU is useless, since they already have the highest result. Higer RTRY values and lower RSTEP values give about 0.5 \% better results at smaller MTU values than 500. For random takeback, seldom takeback and higher RSTEP end up with 2\ldots 4 \% better results at lower MTU than 500 (4 \% for RTRY=12000/RSTEP=5/MTU=400, 2 \% for RTRY=15000/RSTEP=8/MTU=200). But some higher MTU values than 500 also improve the result (3 \% for RTRY=10000/RSTEP=5/MTU=700, 1.5 \% for RTRY=15000/RSTEP=8/MTU=600)

%\subsection{Two-Dimensional Freedom}

%\subsubsection{Fixed takeback at MTU=500}

%In general, the fluctuations are smaller compared to the three-dimensional case (Fig.). Probably the geometric restrictions force frequent resolving of bad moves which have not much time to grow long. Better results are achieved with large RTRY values and also for large RSTEP values. The absolute maximum 290 monomers is at RTRY=15000/RSTEP=10. The control parameters are stable in the ranges of 7000\ldots 16000 for RTRY and 3\ldots 6 for RSTEP. The drawback is that the program efficiency decreases when RTRY increases, thus RTRY=6000/RSTEP=3 can be considered as optimal values (277 monomers, 96 \%).

%%%========== fig:2Dfixtop ==========
%\begin{figure}[!ht]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/peak-unit-p2Dp0205293087r01-f-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:2Dfixtop}
%Relative number of generated monomers with fixed takeback at different RSTEP values and MTU=500 for 2-dimensional freedom}
%\end{figure}

%%%========== fig:2Drandomtop ==========
%\begin{figure}[!hb]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/peak-unit-p2Dp0205293087r01-r-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:2Drandomtop}
%Relative number of generated monomers with random takeback at different RSTEP values and MTU=500 for 2-dimensional freedom}
%\end{figure}
%%%========== fig:2DvariedMTU ==========
%\begin{figure}[h]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/analyys_yle_mtu_2D-RF-unit-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:2DvariedMTU}
%Relative number of generated monomers at selected RTRY/RSTEP values for 2-dimensional freedom}
%\end{figure}

%\subsubsection{Random takeback at MTU=500}

%Clear fluctuations appear upto RSTEP=4 (Fig.~\ref{fig:2Drandomtop}), but they calm down around the maximum result of 290 monomers with parameters RTRY=5000/RSTEP=3. Further increasing of RTRY value at RSTEP=3 does not improve result due to too long resolving time of bad moves. RSTEP=5 shows the most stable behaviour in the range of 4000\ldots 14000 for RTRY, varying less than 5 \% from the local maximum of 274 monomers at RTRY=8000.

%Fixed takeback has 1 \% better maximum result compared to random takeback, still it appears that random takeback resolves more efficiently bad moves without abandoning monomers with good configurations, which may be a case for fixed takeback.

%\subsubsection{Varied MTU}

%The test runs are performed for selected parameters RTRY=6000/RSTEP=3, RTRY=12000/RSTEP=5, \\
%RTRY=14000/RSTEP=4 and RTRY=15000/RSTEP=2 for fixed takeback; RTRY=5000/RSTEP=3, \\
%RTRY=10000/RSTEP=4 and RTRY=14000/RSTEP=6,7 for random takeback (Fig. \ref{fig:2DvariedMTU}).

%A rise of the efficiency by about 8 \% appears at RTRY=15000/RSTEP=2 (fixed takeback) and RTRY=14000/\\RSTEP=7 (random takeback) as a result of MTU decrease from 500 to 300. The result is also 3 \% better when MTU increases from 500 to 700 at RTRY=10000/RSTEP=4 (random takeback). The result for fixed takeback improves at RTRY=12000/RSTEP=5 and RTRY=14000/RSTEP=4 by 2 \% and 1 \% respectively. Rest of the control parameter pairs tested have their maximum at MTU=500.

%\subsection{One-Dimensional Freedom}

%\subsubsection{Fixed takeback at MTU=500}

%The fluctuations due to the small number of monomers in bad moves are even smaller than for the two-dimensional case (Fig.~\ref{fig:1Dfixtop}). Already RSTEP=2 is good enough to use with the absolute maximum at RTRY=9000. At higher RTRY values the efficiency decreases.

%%%========== fig:1Dfixtop ==========
%\begin{figure}[!ht]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/peak-unit-p1Dp0305293087r01-f-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:1Dfixtop}
%Relative number of generated monomers with fixed takeback at different RSTEP values and MTU=500 for 1-dimensional freedom}
%\end{figure}

%%%========== fig:1Drandomtop ==========
%\begin{figure}[!hb]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/peak-unit-p1Dp0305293087r01-r-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:1Drandomtop}
%Relative number of generated monomers with random takeback at different RSTEP values and MTU=500 for 1-dimensional freedom}
%\end{figure}
%%%========== fig:1DvariedMTU ==========
%\begin{figure}[h]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.74\textwidth]{../figures/analyys_yle_mtu_1D-RF-unit-ENG-bw.eps}
%\end{tabular}
%\end{center}
%\caption{\label{fig:1DvariedMTU}
%Relative number of generated monomers at selected RTRY/RSTEP values for 1-dimensional freedom}
%\end{figure}
%\subsubsection{Random takeback at MTU=500}

%RSTEP=2, 3 are quite unstable (Fig.~\ref{fig:1Drandomtop}). There is a plateau in the RTRY range of 12000\ldots 16000 at RSTEP=5 containing absolute maximum of 287 monomers and the efficiency of 89 \% or more. In general, fixed takeback is more stable, but random tackeback is more productive.

%\subsubsection{Varied MTU}

%The selected parameters are RTRY=9000/RSTEP=2,4 and RTRY=13000/RSTEP=3 for fixed takeback; \\
%RTRY=7000/RSTEP=3, RTRY=8000/RSTEP=6, RTRY=9000/RSTEP=4 and RTRY=12000/RSTEP=5 for random takeback (Fig. \ref{fig:1DvariedMTU}).

%Almost all of them have their maximum results at MTU=500. Varying of MTU improves the result only in the case of random takeback at MTU=600 for RTRY=12000/RSTEP=5 and RTRY=7000/RSTEP=3  by 0.5 \% and 1 \% respectively. 

%%=================================
%%========== Conclusions ==========
%%=================================

\section{CONCLUSIONS}

%The optimal values for program control parameters depend a lot on system geometry and applied restrictions. To ensure the desired density in the simulation box with the shortest time some test should be carried on for the specific simulation box geometry as well as for the polymers itself, since the monomer structure and length may have their influence on the generation speed and efficiency; and again, control parameters need some tuning. And one should mention that different initial seed values of the random number generator may speed things up but somevalues do slow down the generation speed. Also, some work is going on to introduce MPI \cite{mpi} into the program to enable the program benefit from the possibilities of parallel computers and to improve calculation speed.

%Since the program is written in modular fashion, it is easily expandable and it is easy to add new features into the input files. Program already outputs data in some common format. The special interest is focused on compatibility with DL\_POLY \cite{bib:dlpoly}, what is a well-known Molecular Dynamics package. Some \textit{mcgen} features are initially inspired from the idea of breaking the generation at the point where the chain is hopelessly trapped and relaxing the chain with Molecular Dynamics simulations, using DL\_POLY. After the relaxation the generation can be continued from the breaking point, but now with different coordinates for atom and ions and therefore with different potential energy of the system.

%%======================================
%%========== Acknowledgements ==========
%%======================================

%%\section*{\Large\textbf
\acknowledgements
%The support and computational resources of the Estonian Grid project are greatly acknowledged.
%%}


%%================================
%%========== References ==========
%%================================

%%\section*{\Large\textbf
%%{References}
%%}
%%\begin{itemize}
\begin{thebibliography}{1}%%{References}

%%\bibitem{baudry}P. Baudry, S. Lascaud, H. Majastre and D. Bloch, 
%%\textit{J. Power Sources}, 1997, \textbf{68}, 432.
%%\bibitem{heitner}K.L. Heitner, \textit{J. Power Sources}, 2000, \textbf{89}, 
%%128.

%%\bibitem{gray}F.M. Gray, \textit{``Polymer Electrolytes''}, The Royal Society 
%%of Chemistry, Cambridge, 1997.

%%\bibitem{croce1}F. Croce, G.B. Appetecchi, L. Persi and B. Scrosati, 
%%\textit{Nature}, 1998, \textbf{394}, 456.

\bibitem{cohen} Y. Bar-Cohen, \textit{Electroactive Polymer (EAP) Actuators as Artificial Muscles - Reality, Potential and Challenges}, SPIE Press, Bellingham, Washington, 2001.

%%%\bibitem{punn01} M. Anton, A. Punning, A. Aabloo, M. Listak, and M. Kruusmaa, ``Towards a biomimetic EAP robot'', in \textit{``Towards Autonomous Systems''}, \textit{Proc. of TAROS 2004}, University of Essex, 2004.

%%%\bibitem{punn02} A. Punning, M. Anton, M. Kruusmaa, and A. Aabloo, ``A biologically inspired ray-like underwater robot with electroactive polymer pectoral fins'', in \textit{Proc. of the Int. IEEE Conf. Mechatronics and Robotics 2004 (MechRob'04)} \textbf{2}, pp. 214-245, Aachen, 2004.

%%%\bibitem{punn03} A. Punning, M. Anton, M. Kruusmaa, and A. Aabloo, ``An engineering approach to reduced power consumption of IPMC (Ion-Polymer Metal Composite) actuators'' in \textit{Proc. of the 12th Int. IEEE Conference on Advanced Robotics (ICAR 2005)}, Seattle, Washington, 2005.

\bibitem{punn04} A. Punning, M. Kruusmaa, and A. Aabloo, ``Surface resistance experiments with IPMC sensors and actuators'', accepted to \textit{Sensors and Actuators: A. Physical}.

\bibitem{THOMAS-SSI-2001-143-83-87} A. Aabloo, and J. Thomas,``Molecular dynamics simulation of lithium ion mobility in a PEO surface'', \textit{Solid State Ionics} \textbf{143}, pp. 83-87, 2001.

\bibitem{THOMAS-JMATCHEM-2003-13-214-218} A. Hektor, M. K. Klintenberg, A. Aabloo, and J. O. Thomas,``Molecular dynamics simulation of the effect of a side chain on the dynamics of the amorphous LiPF$_6$-PEO system'', \textit{J. Mater. Chem.} \textbf{13}, pp. 214-218, 2003.

\bibitem{THOMAS-SSI-2004-168-249-254} H. Kasem{\"a}gi, A. Aabloo, M. K. Klintenberg, and J. O. Thomas, ``Molecular dynamics simulation of the effect of nanoparticle fillers on ion motion in a polymer host'', \textit{Solid State Ionics} \textbf{168}, pp. 249-254, 2004.

\bibitem{THOMAS-CHEMMAT-2005-17-3673-3680} D. Brandell, A. Liivat, A. Aabloo, and J. O. Thomas, ``Conduction mechanisms in crystalline LiPF$_6\cdot$PEO$_6$ doped with SiF$_6^2-$ and SF$_6$'', \textit{Chem. Mater.} \textbf{17}, pp. 3673-3680, 2005.

\bibitem{THOMAS-JMATCHEM-2005-15-1422-1428} D. Brandell, A. Liivat, H. Kasem{\"a}gi, A. Aabloo and J. O. Thomas, ``Molecular dynamics simulation of the LiPF$_6\cdot$PEO$_6$ structure'' \textit{J. Mater. Chem.} \textbf{15}, pp. 1422-1428, 2005.

\bibitem{THOMAS-JMATCHEM-2005-15-4338-4345} D. Brandell, A. Liivat, A. Aabloo, and J. O. Thomas, ``Molecular dynamics simulation of the crystalline short-chain polymer system LiPF$_6\cdot$PEO$_6$ (M$_w\sim$ 1000)'', \textit{J. Mater. Chem.} \textbf{15}, pp. 4338-4345, 2005.

\bibitem{THOMAS-JCOMPHYS-2005-26-716-724}A. Liivat, A. Aabloo, and J. O. Thomas, ``Development of a force field for Li$_2$SiF$_6$'', \textit{Journal of Computational Chemistry} \textbf{26}, pp. 716-724, 2005.

\bibitem{jensen} F. Jensen, \textit{Introduction to Computational Chemistry}, John Wiley \& Sons, Chichester, 2002.

\bibitem{allen}M. P. Allen and D. J. Tildesley, \textit{Computer simulation of liquids}, Oxford University Press, 1986.

\bibitem{rapaport}D. C. Rapaport, \textit{The art of Molecular Dynamics simulation}, Cambridge University Press, 1995.

\bibitem{ercolessi}F. Ercolessi, \textit{A molecular dynamics primer}, Spring College in Computational Physics, ICTP, Trieste, 1997, \textit{http://www.fisica.uniud.it/\~{}ercolessi/md/md/}.

\bibitem{ohno}K. Ohno, K. Esfarjani and Y. Kawazoe, \textit{Computational Materials Science. From ab initio to Monte Carlo methods}, \textit{Springer Series in Solid-State Sciences} \textbf{129}, Springer-Verlag, Berlin, 1999.

\bibitem{taopang}Tao Pang, \textit{An introduction to Computational Physics}, Cambrigde University Press, 2003.

\bibitem{leach}A. R. Leach, \textit{Molecular modelling. Principles and applications}, Addison Wesley Longman Ltd., Reading, Mass., 1998.

\bibitem{stuart}B. Stuart, \textit{Polymer Analysis}, John Wiley \& Sons, Chicester, 2002.

\bibitem{METROPOLIS-JCHEMPHYS-21-1953-1087-1092} N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, ``Equation of state calculations by fast computing machines'', \textit{J. Chem. Phys.} \textbf{21}, pp. 1087-1092, 1953.

\bibitem{colbourn}E. A. Colbourn, \textit{Computer Simulation of Polymers}, Longman Scientific and Technical, Chichester, 1994.

\bibitem{TADOKORO-MACMOL-6-1973-672-675} Y. Takahashi, and H. Tadokoro, ``Structural studies of polyethers, (-(CH$_2$)$_m$-O-)$_n$. X. Crystal structure of poly(ethylene oxide)'', \textit{Macromolecules} \textbf{6}, pp. 672-675, 1973.

\bibitem{bib:gray} F. M. Gray, \textit{Polymer Electrolytes}, The Royal Society of Chemistry, Cambridge, 1997.

\bibitem{WRIGHT-ELCHIMACTA-43-1998-1137-1143} P. V Wright, ``Polyner electrolytes -- the early days'', \textit{Electrochimica Acta} \textbf{43}, pp. 1137-1143, 1998.

\bibitem{ARMAND-SOLSTATION-63-1994-309-319}M. Armand, ``The history of polymer electrolytes'', \textit{Solid State Ionics} \textbf{69}, pp. 309-319, 1994.

\bibitem{THOMAS-JCHEMPHYS-101-1994-10064-10073} S. Neyertz, D. Brown, and J. O. Thomas, ``Molecular dynamics simulation of crystalline poly(ethylene oxide)'', \textit{J. Chem. Phys.} \textbf{101}, pp. 10064-10073, 1994.

\bibitem{mpi} \textit{The Message Passing Interface (MPI) standard}, Mathematics and Computer Science Division at Argonne National Laboratory, Argonne, Illinois, \textit{http://www-unix.mcs.anl.gov/mpi/index.htm}.

\bibitem{bib:dlpoly}
W. Smith, and T. Forester, \textit{The DL\_POLY Project}, TCS Division, Daresbury Laboratory, Daresbury, Warrington WA44AD, England.

\end{thebibliography}

\end{document}
