Molecular Dynamics Simulation of Synthetic Polymers

The polymeric systems are great interesting both academic and industrial level, due to interesting applications in several areas. Nowadays, the trends for the development of new materials are focused on the study of macromolecular complex systems. A macromolecular complex system is formed by the interaction between the macromolecules of the same system. There are many macromolecular complex systems; this chapter will refer to complex systems of synthetic polymer blends(Donald R. P., 2000) and synthetic polymer at the airwater interface(Gargallo et al., 2011). The most synthetic polymers are characterized by flexible structure able to adopt multiple and variable forms. The study about molecular behavior in polymers has attracted the attention of many researchers over the years. Studies in dilute solutions of macromolecules have contributed to solving the problem relating to the molecular characterization. Despite, the emergence of characterization methods and theories that have allowed for interpretation of experimental behavior, researchers still have not fully resolved the problem with interactions at the atomic level in polymer systems. The most polymers are amorphous so in some cases are impossible to obtain information about atomic structure through x-ray techniques or atomic force microscopy. At this point, the computer simulations have played an important role in microscopic understanding of dynamical properties of some polymeric systems.


Introduction
The polymeric systems are great interesting both academic and industrial level, due to interesting applications in several areas. Nowadays, the trends for the development of new materials are focused on the study of macromolecular complex systems. A macromolecular complex system is formed by the interaction between the macromolecules of the same system. There are many macromolecular complex systems; this chapter will refer to complex systems of synthetic polymer blends (Donald R. P., 2000) and synthetic polymer at the airwater interface (Gargallo et al., 2011). The most synthetic polymers are characterized by flexible structure able to adopt multiple and variable forms. The study about molecular behavior in polymers has attracted the attention of many researchers over the years. Studies in dilute solutions of macromolecules have contributed to solving the problem relating to the molecular characterization. Despite, the emergence of characterization methods and theories that have allowed for interpretation of experimental behavior, researchers still have not fully resolved the problem with interactions at the atomic level in polymer systems. The most polymers are amorphous so in some cases are impossible to obtain information about atomic structure through x-ray techniques or atomic force microscopy. At this point, the computer simulations have played an important role in microscopic understanding of dynamical properties of some polymeric systems.
In this chapter a review of molecular simulation methodologies focused on polymer systems will be discussed. This review is based on two case studies. The first one corresponds to a monolayer of amphiphilic polymers and their study at the air-water interface . Three polymers were selected Poly(ethylene oxide) (PEO), Poly(tetrahydrofuran) (PTHF) and Poly(-caprolactone) (PEC) and their cyclodextrin inclusion complexes. The second system corresponds to study the motion of side chains of polymers in condensed phase. For this reason, a diblock copolymers of Polystyrene-block-poly(t-butylacrylate) has been considered (Encinar et al., 2008). The methodology to build the models in both systems will be provided, as well as, the methodology to obtain the atomic parameters, necessary to molecular simulations.
To carry out molecular dynamics studies of polymer systems is necessary to know some important physicochemical properties or thermodynamic parameters of polymers. For the www.intechopen.com other hand, it is important to know the different methodologies employed to carry out molecular dynamic simulation of polymers.

Physicochemical properties of polymers
In this chapter we concéntrate efforts on giving the reader a brief description of some important physicochemical properties of the polymers for molecular simulation analysis analysis. Macromolecular conformation is an important point to consider. The macromolecular chains in dilute solution have variables internal rotational angles, Figure  1. The internal rotational angle  can take values that fluctuate within a range so that different minimum energy can be obtained for various values of the angles . The energy difference associated with the variation of this angle is minimal. This leads to a large number of stable conformations and significant mobility to oscillate around these conformations. Both factors cause that the macromolecule has a significant flexibility. Experimentally the dimensional parameters that characterize a macromolecule can only be characterized in terms of statistical parameters that represents an average over all conformations. The end-to-end distance and radius of gyration are statistical parameters that characterize the spatial distribution of a chain macromolecule (Flory, 1969(Flory, , 1995. The end-to-end distance represents the distance between end chains. This parameter can vary from a maximum value, where the polymer chain is fully extended to a minimum value corresponding to the sum of the van der Waals radii of the terminal groups, Figure 2. The end-to-end distance is a statistical parameter that characterizes the spatial distribution of the chains relating the microscopic with macroscopic properties, also allows to relate the hydrodynamic parameters with thermodynamic parameter, and finally allows to obtain information respect to flexibility chain.  (1) where W(r) is probability distribution function, defined by But the chains may adopt different conformational states depending on the characteristics of the chain, features such as substituents and the solvent. For the different conformational states, there are different chains models. Is not intended to provide a thorough review of statistical mechanics of polymer chains (Flory, 1969), but is important to note that the value of depends on two main factors, "short-range" and "long-range" factors. For these two factors, different chain models were raised in the past. Short-range factors involved characteristics of the bonds and interaction between neighboring atoms or groups along a chain, this chain correspond to a rotation free model. Hence, can be calculated knowing the values of bond angles and bond distances for all the conformations, the value obtained is called dimension mean square of end-to-end distance . Take into account factors such as bond angle restrictions and steric hindrance between substituent groups, we can obtain unperturbed dimension , and may be expressed by equation 2: (2) where  is the stiffness coefficient or conformational factor. The is so-called "theta condition".
In long-range factors are involved the interactions between distant segments along the chain, as well as interaction with the solvent. Taking into account that, two distant segments of a chain, cannot occupy the same place at the same time, arise the excluded volume concept. Under effect of these interactions the end-to-end distance to linear chains in dilute solution, may be described by: where  linear expansion chain factor. In equation 3, the root-mean-square end-to-end distance depends of number of segments, temperature and solvent. When >1 them the interaction between solvent and chains are favorable, this means that excluded volume effect is important. Instead, if the interactions between segments and solvent are unfavorable them  tend to decrease; this causes folding of the chains to minimize the contact with the solvent. These effects can be balance and them the polymer chain adopts an unperturbed state.
The radius of gyration is the parameter more appropriate to describe the dimensions of a macromolecule. The mean-square radius of gyration represents the mean square of the distance between each atom in the chain and center of mass of the chain and may be defined to following equation 4: where n is the segments numbers of the chain and represents the mean-square distance of the i-th segment from the center of mass of the chains. Through the random flight model we can express the radius of gyration in terms to mean-square end-to-end distance, indicating their interchangeability. For a freely rotating chain of n segments of length l, the following relationship is established: A complete derivation of these equations can be reviewed in various texts of statistical mechanics of polymers (Flory, 1995) (Yamakawa, 2001). We have said that the dimensional parameters are important to know because of their contribution to the characterization of polymers. In addition, these parameters can be calculated both experimentally and theoretically. For example, radius of gyration has been calculated, in several works, from molecular simulation trajectories and compared with the experimental values (Lee et al., 2006;Prosa et al., 2001Prosa et al., , 1997.

Molecular dynamic simulation of supramolecular structures at the air-water interface
A great interesting in the last year toward to study theoretical of self-assembly of macromolecular systems have been report (Duncan et al., 2011;Korchowiec et al., 2011;Saldias et al., 2009). The self-assembly of macromolecular systems are great interesting due to applications in several areas, as medicine, drug delivery and smart materials (Wick and Cummings, 2011). The polymers at the air-water interface are considering self-assembly systems. The studies of polymers at the air-water interface considering several molecular simulation methodologies, such as atomistic ) and coarsegrained (Duncan et al., 2011) models. In this chapter we will focus on atomistic molecular simulation methodology. The aim of studies of polymer at the air-water interface, using molecular simulation, is characterized the orientation of polymer chains at the interface, as well as, calculate the different interaction between amphiphilic polymer and the water subphase. www.intechopen.com

Build of the molecular models
In general, the models used to reproduce several polymer chains at air-water interface, consist of the build a bilayer of polymers with an enough thick water slab between of the chains. All the systems were simulated taking into account three experimental surface areas 10, 20 and 30Å 2 per molecules. The water slab should be larger in one axis (i.e. z-axis), to avoid the contact between the polymer chains, see Figure 3. A vacuum space larger enough to avoid the contact between chains of different side of the bilayer, should be considered. To the constructor of different molecular structure (repeating unit of the polymers) any molecular builder is required, the important point is generated the pdb coordinate file. To generate the polymer chains a script was made in Tcl language and them executed in VMD (Humphrey et al., 1996) program and this way, the polymer chain of ten or twelve repeating unit were obtained. The build of water slab was made to using VMD program. For this purpose, was necessary obtained the coordinate of each polymer chain and them located these chains in the correct coordinates in order to build the water slap in the center of the coordinate and this way, obtained the model that represent several polymer chains oriented at the air-water interface, as  To the build of the inclusion complexes of cyclodextrin-polymers,  and -cyclodextrin were considering. The cyclodextrin structures were obtained from Protein Data Bank, PDB entry 1CXF for -CD and 1D3C for -CD. These structure were fully minimized, using CHARMM force field and a good agreement with the experimental geometry was obtained. To carry out the minimization and molecular dynamic simulation of synthetic polymer was necessary to use a force field, which content the atomic parameter appropriate to a particular polymer. To obtained the atomic parameters to Poly(ethylene oxide) (PEO), Poly(tetrahydrofuran) (PTHF) and Poly(-caprolactone) (PEC) we used the CHARMM force field, this force field is open source and it is possible aggregate the parameter because the polymers implicate in the study was necessary calculate the atomic parameter, through quantum mechanics. For this purpose, Paratool package implemented in VMD (Humphrey et al., 1996) software was used. Six models were built, three corresponding to the inclusion complexes: PEC-CD, PEO-CD and PTHF-CD and the polymer respectively at the air water interface, see Figure 4.

Considerations to molecular dynamic simulations
All Molecular dynamic (MD) simulations were carried out using NAMD (Kalé et al., 1999) program and the CHARMM27 (Kuttel et al., 2002)   NPT ensemble; the non-bonded van der Waals and electrostatic interactions were calculated with a cut-off using a switching function starting at a distance of 10Å and reaching zero at 12Å. The TIP3P (Jorgensen et al., 1983) water model was employed for the solvent. Periodic boundary conditions were applied with a rigid cell. The particle mesh Ewald (Darden et al., 1993) (PME) method was employed for computation of electrostatic forces. An integration time step of 1fs was assumed, permitting a multiple time stepping algorithm (Grubmüller et al., 1991) to be employed in which interactions involving covalent bonds were updated every time step. In all simulations, Langevin dynamics were utilized to keep a constant temperature of 298K, likewise the hybrid Nosé-Hoover Langevin piston method (Izrailev et al., 1997) was used to control a constant pressure of 1 atm. The system was first equilibrated for 1ns. The backbone carbons of the polymers in the inclusion complexes were restrained with a 0.5kcal/mol/Å 2 spring constant during relaxation. Then, simulations of 5ns were performed to all the systems. A novel alternative to obtain the charges and potentials to each atoms or groups is through CHARMM General Force Field (CGenFF) (Vanommeslaeghe et al., 2010). This force field is an extension of the CHARMM force field to small molecules and included a wide range of chemical groups. The molecule, in our case the monomer, is required in mol2 format to be uploading in the web site of Paramchem and this way, obtain the charges and potentials atomic. The server provides the topology and parameter files, necessaries to run the MD simulation. Topology file contains the monomer residue, atom types and the charges. Then, is mandatory to make changes to the topology file because we need create a topology for a polymer containing many repeat units. At the moment the topology file only have one residue that corresponding to the monomeric structure, so is require to build three types of residues. A "head" residue that defined the initial monomer chain, the "repeat" residue that is corresponding to the residue in-between chain and "tail" residue that defined the final monomer chain. Another way is to define only one residue corresponding to de monomer and build a path union, where is specified the hydrogen atoms that most be eliminated and the atoms that most be linked. This way is possible obtaining the coordinate and connectivity file to a polymeric chain.

Properties polymeric systems, calculations and analysis
From MD simulations trajectories several interesting properties cab be obtained and give insight about the behavior of inclusion complexes and polymer chains at the air-water www.intechopen.com interface. The most important point to know about the supramolecular systems (Komarov et al., 2011;Krishnamoorthy et al., 2011) at the interface is the orientational correlation about the order of monolayer. In the case of study of the inclusion complexes (Almeida et al., 2011;Boonyarattanakalin et al., 2011;Pan et al., 2011;Shi et al., 2011), we calculate the order parameter (Sr). This property is defined in terms of second-order Legendre polynomials and is the average of an angle made by a section of the polymer chain in a specified direction. This orientation correlation function S(r) is defined by equation: .5 3 where  is the angle between i and j vectors, when the center of mass of the associated chains of the vectors i and j are separated by the distance r. For the case of the inclusion complexes, we selected the atoms carbon C-C unit vectors of the different chains from PTHF, PEC and PEO as the vectors i and j. This parameter is indicative of the tendency towards polymers ordering in systems containing flexible chains. A value of 1.0 for S(r) mean that the chains are parallel to each other and cero value for a disordered distribution.
In the Figure 5, the S(r) variation as a function of time for these polymers is shown. It is clear that PEC is more organized at the air-water interface than PTHF and PEO. These results were in good agreement with the experimental data. Comparing the order parameter for PEC, PTHF and PEO and their respective ICs, can observe that the ICs are more ordered at air-water interface than their precursor polymer, Figure 6 a, b and c.
The order parameter shown differences evident in lose conformational freedom for cyclodextrins-polymers systems. The inclusion complexes are stiffer than single chains explaining the lower values of S(r) of the chains compared to those included in cyclodextrins.
The radial distribution function, g(r), is another important property that describe the probability to find an atom in a distance r from another atom and can be represent by equation 7 Figure 7 and 8 presents g(r) between the hydrophobic groups of the different polymer chains and also g(r) between hydrophilic groups of the different polymer chains and water molecules. The electrostatic interactions between hydrophilic groups of the precursor polymers and water molecules are established. Nevertheless, some differences were seen in g(r) plots to PEC, PEO and PTHF. In the case of, PEC and PTHF (Figure 7a,b) the distance between the hydrophilic groups and water molecules are slightly lower in compare to the distances in PEC.
The parameters can be estimate from molecular dynamics simulation trajectories, in the same way that the radial distribution function and order parameter. Scripts for each of the parameter were written in Tcl language and were ran in VMD program. The scripts can be to provide the reader through a mail to the author. www.intechopen.com

Free energy conformational calculation to the copolymer systems at condensed phase
The di-block copolymers in melt state have been studied by divers researcher groups, due to great interesting from the point view of the dielectric and dynamic-mechanics behavior (Alig et al., 1992;Encinar et al., 2008;Karatasos et al., 1996;Kotaka and Adachi, 1997;Kyritsis et al., 2000;Ma et al., 2003;Moreno and Rubio, 2002;Vogt et al., 1992). Dynamic-mechanical spectroscopy is an important technique to characterize the deformation of polymeric materials applying, for example, a tension to material. The viscoelastic behavior of the polymers can be studied through the determination of properties such as the complex modulus and glass transition temperature. On the other hand, the dielectric relaxation explore the response of a sample submitted an electric field. The permanent dipolar moments of a flexible polymer can be reoriented by an electric field, giving the relaxation mechanics of the fluctuations of dipolar moieties of the sample. For this reason is interesting to know at atomic level the conformational motions that have a polymer or copolymer. The conformational free energy difference (G) for a polymer chain, can be calculated through molecular dynamics simulation at different temperatures, with the aims of obtain a conformational wide sampling. For this purpose, Boltzmann equation (Gilson et al., 1997) was used to obtain G, equation 8.

∆ (8)
Where f is the frequency distribution of the conformational change of the dihedral angle of side chain of PtBa, R is the gas constant and T is the temperature.
The case study was carried out for two systems of diblock copolymer of Polystyrene (PS) and Poly(t-butylacrylate) (PtBa) with different Poly(t-butylacrylate) chain lengths. The two blocks of PS and PtBa are amorphous and immiscible, also these blocks have well separates their glass transitions temperature (Tg). The PS-b-PtBa copolymer is dielectrically active due to PtBa part only. Nevertheless, the PtBa block has no present dipole moment to along the chain, thus the segmental mode is not influenced by other dynamical contributions. Besides, both blocks are active in dynamic-mechanical spectroscopy. The results obtained from the two techniques were compared with the PtBa homopolymer. The experimental results shown that the shape of the relaxation curves has different dependences with the temperature in the copolymers and the PtBa homopolymer. However, the relaxation maps to -mode were almost equivalent in the copolymer than in the PtBa homopolymer, but the relaxation time to -mode were different in the copolymer respect to homopolymer. This result was unexpected, because the -relaxation has a local character due to the PtBa side chains relaxation. The above behavior could be owing to the influence of the PS on the PtBa side chains. For this reason, different environments were simulated in condensed phase and three models were built Figure  Paratool is a tool that used quantum mechanics methodology to the optimized structures and vibrational frequency to obtain the atomic parameter, both unbounded and bounded generating topology and parameter files, in CHARMM format. The MD simulations were carried out in NAMD program using CHARMM force field. The polymeric fragments containing five repeating units of each sequence were placed in a box under periodic boundary condition Figure 10; ten chains were placed in three different boxes for each model. Molecular simulation of 1ns of long time were carried out using Valbond Force Field (CVFF) (Maple et al., 1998) and CHARMM forcé field. The three models were run at different temperatures between 1000K and 1700K with the purpose of obtain a wide conformational sample. A conformational analysis was performed to each model in order to estimate the free energy difference between two conformations for dihedral angles  1 and  2 , see Figure 9.
The frequency average of the dihedral angles was calculate in the three models at the different temperatures from molecular simulation trajectories Figure 11, and applied to the equation 8. The G values were compared with activation energy E A from dielectric experiments to the -relaxation, Table 1.   The experimental data showed that the characteristic time to -transition is faster in the case of homopolymer compared to the copolymers. This behavior provides evidence that the PS blocks are influencing to the PtBa blocks. From MD simulations was possible to obtain the frequencies distribution to the conformational changes of 1 and 2 dihedral angles, Figure 11. Distribution frequencies calculated for P(305), P(307) and PtBa for all temperatures were applied to equation 8 and them the G was obtain from the plots of ln(f 1 /f 2 ) against 1/T. Conformational free energy had the following trend PtBa>P(305)>P(307). Table 1 shown that the G is unfavorable in the case to PtBa, indeed this mean that there is less mobility of ester group, then the relaxation processes should be more hindered than in the other fragments. These results explain the activation energies values obtain by experimental techniques, where the E A to PtBa is higher in compare to P(305) and P(307). Therefore, the results obtain by MD simulation of different fragment of the copolymer are in agree with the experimental data. In addition, the model containing the monomer surrounded by PS and PtBa (P305) have a higher mobility, so the activation energy to the relaxation process is less.
We have seen that the molecular simulation methods can be effective to study behavior of the system in solution as well as in condensed phase. The considerations to carried out molecular simulation in polymeric systems are: 1. To know if the process of study occur in solution or condensed phase 2. The time scale that take this process, so with that information is possible decide which methodology can be better to establish an equilibrium between simulation time and accuracy of results. Some process, as the mesoscopic structures, should be studied with coarse-grained methodology. 3. Establish that force field is more appropriate to our polymeric system. 4. In the absence of atomic parameter is necessary get through quantum mechanics and add to the force field, in the case of those of open source, as CHARMM. 5. The atomic parameter can be to get of an easy way through Paratool module available in VMD program or using the web platform of Paramchem. 6. Building models should represent the issue to be solved. 7. Select to appropriate program to carried out the MD simulation 8. Depending of the issue of our system, may be necessary apply to the polymer system constrain. Oftentimes, is required the study of movement of the side chain and not necessary to considered the backbone polymeric movement, in this case can be considerer apply restriction to the movement of the backbone. The restraints to the motion have an important advantage that allows simulation time economy. The important is emphasizing de differences between restraints and constraints. The constraints eliminate the hard degrees of freedom of the systems that correspond to high frequencies of vibration. The constraints can be applied on the hydrogen bond distances and hydrogen bond angles, for this purpose the shake algorithm is used.

www.intechopen.com
Shake algorithm (Andersen, 1983;Ryckaert et al., 1977) to do cyclic computation of all constraints, in an iterative fashion until each constraint be satisfied. 9. Other input parameter to considerer before starting with de MD simulation is the canonical ensemble. Simulation to the polymer at the air-water interface the NPT ensemble is appropriate to take into account. 10. Finally is important the simulation time, it is necessary that the polymer chains to reach equilibrium. For this reason, the molecular simulation energy should be plotted against the time simulation; if the energy is not constant it will require longer simulation.

Conclusions
Molecular dynamic simulations can became a powerful tool to study of the behavior of polymer system both in solution and condensed phase. In this chapter, two different works were approached, in both cases the intermolecular interactions were calculated for the purpose of established a correlation between theoretical and experimental results. In the case of inclusion complexes, molecular simulation was carried out to understand the behavior at the air-water interface. The results obtained from molecular simulation trajectories to each system were used to calculate the order parameter and correlated with -A isotherms from Langmuir experiments. Order parameter established that both PEC and IC-CD-PEC were more ordered systems related to PTHF and IC-CD-PTHF and PEO and IC-CD-PEO. These results were agreement with experimental results. However, radial distribution function was calculated with the purpose of obtain information about orientation of polymeric systems at the interface. From radial distribution function could be established that the interaction between polar groups of the polymers and water molecules took place around 2Å distance, then was possible estimate that there is a diffusion of the polymer to the water. Also, radial distribution function between hydrophobic groups of the different chains was calculated. The distances between polymeric chains were considerably short so the chains adopted a conformation rather widespread; with the purpose of interact among them. In the second case of study, molecular dynamics of Polystyrene-bock-Poly(tbutylacrylate) (PS-b-PtBa) copolymers were carried out at different Polystyrene composition for the purpose of observed the influence of PS in the experimental dielectric relaxation of PtBa. From molecular simulation trajectories the conformational free energy was calculated to the different PS composition in the copolymers. Conformational free energies calculate from Boltzmann equation were correlated with activation free energies calculated from experimental dielectric relaxation. Activation free energy to PtBa was higher compared to P(305) and P(307), this mean that a more hindered change conformational was observed. A less favored change conformational was obtained to PtBa compared to P(305). Therefore, good agree was obtained between the activation free energy results and conformational free energy from molecular simulation. Finally, it is possible conclude that molecular dynamics simulation is a important and challenging tool that can be used to know properties and processes that occurs in polymer systems both dissolution and condensed phase. A important point to take in account is dispose of appropriate atomic parameters, that mean an adequate force field and the other hand, a model that represent the polymeric system.