MCDP: An Advanced Tool to SimulateComb-Like Polymers

SALVADOR LEÓN,1 CARLOS ALEMÁN,1 FRANCESC ESCALÉ,2

MANUEL LASO3

1Departament d’Enginyeria Química, ETSEIB, Universidad Politécnica de Catalunya, Diagonal 647,Barcelona E-08028, Spain2Centre Europeu de Paralellisme de Barcelona, Universitat Politécnica de Catalunya, Jordi Girona 1-3,Barcelona E-08034, Spain3Departamento de Ingeniería Química, Universidad Politécnica de Madrid, José Gutierrez Abascal 2,Madrid E-28006, Spain

Received 24 March 2000; accepted 14 July 2000

ABSTRACT: A strategy to study polymeric systems with ordered structures,and in particular comb-like polymers, is presented. These are dense systems forwhich atomistic simulations with conventional methods are difficult or evenimpracticable. The strategy, which has been incorporated into a computerprogram named MCDP, is based on a Configuration Bias Monte Carlo algorithmand a method to investigate the structure of crystalline polymers using force-fieldcalculations. To obtain a maximum efficiency, the MCDP computer program hasbeen optimized and parallelized. The ability of MCDP to investigate orderedpolymers have been tested by simulating two complex systems: (1) the crystalstructure of poly(4-methyl-1-pentene), and (2) the biphasic structure ofpoly(α-octyl-β-L-aspartate), a comb-like polyamide derived from β-aminoacids. The results obtained from MCDP simulations demonstrates the efficiencyand reliability of this method to study both the NVT and NPT behavior ofordered dense polymers. c© 2000 John Wiley & Sons, Inc. J Comput Chem 22:162–171, 2001

Keywords: comb-like polymers; ordered structures; Configuration BiasMonte Carlo algorithm; MCDP

Correspondence to: C. Alemán; e-mail: aleman@eq.upc.esContract/grant sponsor: DGICYT; contract/grant number:

PB96-0490Contract/grant sponsor: Comisionat per a Universitats i Re-

cerca de la Generalitat de Catalunya; contract/grant number:1999SGR-00142

Journal of Computational Chemistry, Vol. 22, No. 2, 162–171 (2001)c© 2000 John Wiley & Sons, Inc.

MCDP

Introduction

C omputational methods based on force-fieldcalculations have been an important tool in

the study of the structure and properties of syn-thetic polymers. A particularly interesting problemwithin this field is the simulation of multichaindense polymer systems. Thus, the configurationalsampling of a dense system is quite difficult be-cause the motion of individual entities is hindered.Conventional methods, like molecular dynamicsor simple Monte Carlo (MC) algorithms, exhibitthus a poor efficiency.1 Different approaches havebeen designed to overcome such problem. For in-stance, progress has been made trying to reducethe characteristic time scale of the system by us-ing coarse-graining methods.2 In the last 10 yearsa number of novel algorithms for both lattice andoff-lattice MC simulations of condensed polymershas been proposed.1 – 10 However, in all cases thesemethods were designed to study disordered con-densed phases, i.e., mainly melts and amorphousstructures.

On the other hand, comb-like polymers havebeen subject of recent interest in materials sciencebecause they may form liquid crystals. In particular,polymers constituted by either α- or β-amino acidswith attached flexible side chains have receiveda special attention.11 – 18 The most representativemembers of these families of comb-like polymersare the poly(γ -alkyl-α-L-glutamate)s11 – 15 and thepoly(α-alkyl-β-L-aspartate)s.16– 18 These consist ofpeptide α-like-helices, which impart rod-like char-acter to the macromolecules, coupled with longaliphatic side chains emerging from the repeat units.The helices are arranged in layers, whereas the sidechains are in a separated phase, which is locatedbetween the layers (Fig. 1). These ordered biphasicstructures have been well characterized by different

FIGURE 1. Typical packing structure of a comb-likepolymer.

experimental techniques, but their atomistic simula-tion is extremely difficult.

In this article, we present an efficient strategy forstudying the structure and properties of dense poly-mers with ordered structures, in particular comb-like systems, using off-lattice MC simulations. Thismethodology has been incorporated into a com-puter program named MCDP (Monte Carlo simu-lations of Dense Polymers). The article is organizedas follows. In the next section we discuss the basicelements of MCDP and the program structure. Sub-sequently, the suitability of the program to studydense polymers is presented and critically evalu-ated by performing some test simulations. The lastsection contains the conclusions and a discussion ofpossible applications of the program.

Methods and Results

BASIC ELEMENTS OF THE MCDP PROGRAM

The MCDP strategy requires two basic elements.These are: (1) an efficient off-lattice MC algorithm;and (2) a computational method able to generatethe atomic coordinates of ordered structures andevaluate their energies. The MC algorithm selectedfor this work is the Configurational Bias3, 4 (CBMC),whereas a strategy recently developed by us19 andimplemented in the PCSP (Prediction of the CrystalStructure of Polymers) computer program has beenused to deal with the structures.

Configurational-Bias Monte Carlo

This method was implemented for the off-latticetreatment of realistic polymer systems at the begin-ning of the last decade.3, 4 The CBMC was describedin detail in these references. This algorithm consistsof the following three steps: (1) a chain is selectedat random; (2) the chain is cut at a random position;and (3) the chain is sequentially regrown. The chainis regrown bond by bond by examining a number ofpossible positions (Ns) which can be chosen eitherrandomly, or in a biased way, depending on the tor-sional potential associated to the bond. The weightof the regrown part of the chain corresponds to theproduct of the weights of the regrown bonds, whichin turn, are proportional to its normalized Boltz-mann factor. To ensure the condition of microscopicreversibility, an appropriate transition probability ischosen.1, 3, 4

The CBMC method has been successfully usedto simulate melts of dense systems composed by

JOURNAL OF COMPUTATIONAL CHEMISTRY 163

LEÓN ET AL.

chains of medium length.7 It is a very powerfulmethod when the system contains a large concen-tration of chain ends. This is the case of comb-likepolymers in which a large number of alkyl chainsare attached to the main chain. This method has alsoproven to be efficient for the simulation of tetheredchains, which is the situation encountered in comb-like polymers.20 However, it should be emphasizedthat in the work of Siepmann and McDonald alkylthiol chains were attached to an atomically flat sur-face, while in comb-like polymers the chains areattached to a core with a shape defined by the mainchain conformation. It should be mentioned thatconventional CBMC is not the most appropriatedmethod to simulate systems with branched archi-tectures, alternative methods being proposed in thelast years.21 However, the problems associated toCBMC have been avoided in this work by usingfixed bond lengths and angles.

The PCSP Strategy

This is a computational method recently devel-oped to predict and analyze the packing in crys-talline polymers.19 The strategy is based on empir-ical force-field calculations, and evaluates all themodels of packing for a given lattice dimensions,being able to predict the lowest energy one. Forthis purpose, a crystal lattice is generated using theunit cell parameters, which are usually obtainedfrom X-ray and/or electron diffraction, and both thenumber and position of chemical repeating units(CRUs) in the unit cell. Then, the different pack-ing arrangements are obtained by varying the dis-placement between different sheets along the a-, b-,and/or c-axis, as well as the setting angle betweendifferent chains. The energy of each arrangement iscomputed by adding the bonding and nonbondingcontributions, i.e., torsional, electrostatic, and vander Waals contributions, of the chains contained inthe unit cell. As a result of the lattice periodicity,the nonbonding interactions are computed not onlywithin the chains of the unit cell, but also betweena reference unit cell and external units. This methodhas demonstrated to give very good results for bothconventional and modified nylons.19, 22 – 24

Despite its power, the PCSP strategy has an im-portant shortcoming, namely, chain conformationis not relaxed for every packing arrangement. Theconsideration of fixed chain conformations has im-portant advantages from a computational point ofview. However, it reduces the applicability of themethod to compounds with rigid conformationsand without flexible side chains.

PROGRAM STRUCTURE AND IMPLEMENTATION

The general flow chart of the MCDP program isdisplayed in Figure 2. MCDP uses as INPUTs: (1) thegeometry and connectivity of the CRU, (2) struc-tural parameters related to the lattice dimensions,the helical symmetry of the chain, and the positionof the chains within the unit cell; (3) the force-fieldparameters, and (4) user information and direc-tives about the MC simulation. In what follows, amolecular chain is generated and oriented with re-spect to a Cartesian coordinate system such that thechain axis is along the c-axis. Then, the crystal lat-tice is built using the information provided by theuser. The atomic coordinates of the first chain areused to generate all the other chains by applyingthe required translation and/or rotation transfor-mations.

The program allows variance of the following pa-rameters: (1) the conformation of the main chain,(2) the conformation of the side chain, (3) the cell pa-rameters, (4) the setting angle, and (5) the positionsof the chains within the cell. The types of parame-ters to be moved as well as the fraction of movesof each type are chosen by the user according tothe requirements of the simulation. The Metropo-lis algorithm,25, 26 which is the simplest MC scheme,is suitable for all these parameters with exception

FIGURE 2. General flow chart of the MCDP program.

164 VOL. 22, NO. 2

MCDP

of (2). Thus, these parameters, in particular (3), (4),and (5), are usually determined by X-ray and/orelectron diffraction, and therefore, only small oscil-lations around the experimental values are desir-able. There are a number of computational strate-gies able to determine the conformation of the mainchain in terms of atomic coordinates taking ad-vantage of the experimental information.27 – 30 Thisconformation, which is usually assumed to be quiterigid in ordered structures, should be used as start-ing point in MCDP simulations.

On the other hand, the aim of this study is the de-sign of a computational tool for studying comb-likepolymers. Accordingly, the configurational sam-pling of side chains should be much faster than thatof the main chain. This is achieved by using theMetropolis algorithm to sample the degrees of free-dom of the main chain. It should be emphasizedthat, despite the main chain moves, the helix sym-metry can be optionally retained to the experimen-tally observed value along all the simulation. Thisis achieved by the following algorithm: (1) a mainchain torsional angle is selected at random; (2) theMetropolis algorithm is applied to moved the re-maining torsional angles; (3) the value of the se-lected torsional angle is chosen to satisfy the helicalparameters imposed by the helix symmetry.27 It isworth noting that when more than one value of theselected torsional angle satisfy such parameters, thevalue closest to the original one, i.e., that of (1), ischosen. Finally, the CBMC algorithm is used to ex-plore the configurational space of the side chains,although a fraction of Metropolis moves can be alsoused. Figure 3 shows the flow chart of the subrou-tine used in the program to carry out the CBMCmoves of the side chain.

The energy of the system is estimated by consid-ering the torsional and nonbonding contributions.The latter is obtained by adding the van der Waalsand electrostatic terms. The van der Waals energy iscomputed in the usual pairwise additive way usinga Lennard–Jones 6–12 potential, which is only cal-culated between atoms of different chains or atomsin the same chain separated by at least three bonds.The electrostatic term is computed by applying theCoulombic expression. All the nonbonding inter-actions within atoms of the same chain separatedby exactly three bonds (one to four interactions)are reduced by applying of a 0.5 scale factor.31 Thetorsional energy is computed using the classicalthree-term Fourier expression. Bond lengths and an-gles are kept fixed at their equilibrium values alongthe simulations.

FIGURE 3. Flow chart of the subroutine used in theprogram for the CBMC moves of the side chain.

PARALLELIZATION AND SCALABILITY

Parallelization was performed on an Origin 2000Silicon Graphics with 64 MIPS processors at the“Centre Europeu de Paralelisme de Barcelona”(CEPBA). Before parallelizing the code, the first stepwas to find the most computation intensive partof the MCDP program. An analysis of the profilesfor different runs indicated that the routine whichperforms the CBMC moves takes more than 85%of the execution time. Thus, we decided to par-allelize such routine by modifiying the sequentialversion of the code and using the SGI OpenMPcompiler. A number of trials were performed con-sidering both static and dynamic schedulings. Theuse of static scheduling provided a low speedupbecause the work was not well balanced betweenthe different threads. Conversely, the work wasrightly balanced using dynamic scheduling. How-ever, the dynamic scheduling spends some time tosyncronize the threads.

The best speedup was obtained using staticscheduling with chunks of two iterations. Thisscheduling, called interleave, allows one to obtainthe work well balanced because it drecreases in

JOURNAL OF COMPUTATIONAL CHEMISTRY 165

LEÓN ET AL.

FIGURE 4. Profile obtained from a run of the parallelized version of the MCDP program. The parallel regions arerepresented in gray. Note that the work is rightly balanced between the threads. The sequential regions are representedin white (starting and ending parts).

linear form between iterations. The parallelizationof the code using this procedure is illustrated inFigure 4. As can be seen, all the work runs in par-allel (gray color) with exception of the starting andending parts (white color), which run in sequence.Indeed, the sequential parts are so short that it isdifficult to detect them in Figure 4.

The scalability of the MCDP program was inves-tigated by performing some benchmarks.The resultsfor a work of 500 CBMC steps are displayed in Ta-ble I. It is worth noting that very good speedupswere obtained up to 16 processors. For more than16 processors the speedups were not able to growwell.

TEST CASES

The ultimate goal of the MCDP strategy is to pre-dict the structure of ordered polymeric systems, inparticular of comb-like polymers. To investigate thereliability of the MCDP program, test calculationswere performed on some prototypical systems.

TABLE I.Performance of the MCDP Program.

Version No. of Processors Time Speedup

Sequential 1 346.70 —Parallelized 2 177.41 1.95Parallelized 4 89.68 3.86Parallelized 8 45.48 7.62Parallelized 16 23.46 14.77Parallelized 24 17.40 19.92Parallelized 32 14.25 24.33Parallelized 48 11.24 30.84Parallelized 60 10.20 33.99

Time required (in seconds) by the different versions of theprogram to run 500 steps of CBMC of a typical comb-likepolymer. The time and the speedup of the parallelized versionis given as a function of the number of processors. Calcula-tions were run in a SGI Origin 2000 server.

NVT Simulations of CrystallinePoly(4-methyl-1-pentene)

Isotactic poly(4-methyl-1-pentene), abbreviatedPMP, is a semicrystalline polymer with rather com-plicated packing of the bulky pendant group.32 Itscrystal structure consists of four 7/2 helices packedin a tetragonal unit cell of dimensions a = 18.70 Åand c = 13.68 Å. This structure was investigated atatomic level by refining the X-ray data with a con-strained least-squares (CLS) method. Unfortunately,the structure was not successfully determined, dueto both the poor diffraction data and the rather com-plicated structure, i.e., four molecules packed in alarge unit cell.

The performance of the MCDP strategy to studypolymers with bulky side chains has been testedby investigating the packing of PMP. The startinggeometry for the MCDP simulations was taken fromRef. 32. Accordingly, four helices were packed in atetragonal box at the positions indicated in Figure 5.In the crystal, it was found that at each position thehelix can be upward (U) and downward (D) withprobability 1/2. It was considered that the arrange-ment displayed in Figure 5 is the most suitablediscrete distribution of a statistical disorder with re-spect to the U and D sense of the four 7/2 helices.

Two different NVT simulations were performed.In the first, that will be denoted MC/1, the torsionalangles of the helix backbone were kept constant atthe values reported in Ref. 32, the degrees of free-dom being the torsional angles of the alkyl sidechain and the setting angle of each helix. In thesecond simulation, denoted MC/2, the degrees of

166 VOL. 22, NO. 2

MCDP

FIGURE 5. Representative microstructure of the PMPobtained in the MC/1 simulation. U and D refer to theupward and downward sense of the 7/2 helices,respectively.

freedom were those of MC/1 and, additionally, thetorsional angles of the backbone. In all cases the dis-tance between helices within the box was fixed atthe experimental values. The simulations were per-formed at T = 298 K considering both periodiccontinuation conditions and the minimum-imageconvention. Residue-based cutoff were applied at10 Å for the nonbonded interactions. The parame-ters for the different atoms were taken from theAmber force field.31

One of the microstates obtained from the MC/2simulation is displayed in Figure 5. On the otherhand, Figure 6 shows the energy of the systemagainst the number of steps for both MC/1 and

FIGURE 6. Energy of the system as a function of thenumber of steps for the MC/1 and MC/2 simulationsof PMP.

MC/2 simulations. It is worth noting that for thetwo simulations the energy of the system relaxesrapidly. Another interesting feature can be observedin Figure 6: the system becomes much more sta-ble in the MC/2 simulation than in the MC/1 one.Thus, an energy gain of about 250 kcal/mol isachieved when the torsional angles of the back-bone are allowed to move. A detailed analysis ofthe microstructures generated from the two simula-tions indicates that the stabilization of the packing ismainly achieved by the relaxation of the backbone.

Figure 7 shows the distribution of torsional an-gles obtained from MC/1 and MC/2 simulationsfor the seven residues contained in the repeat ofone helix of PMP. The population analysis has beenperformed considering the two torsional angles ofeach side chain as a pair. Thus, the conformationshave been grouped in three sequences: trans-trans,trans-gauche+, and the remaining conformations. Itis worth noting that the results obtained from thetwo simulations are almost identical. Thus, thetrans-gauche+ is the predominant sequence with apopulation of about 98%. The same results havebeen obtained for the remainder three helices con-tained in the unit cell. It should be mentioned thatthe trans-gauche+ sequence was also the conforma-tion achieved by CLS refinements.

On the other hand, Table II compares thetorsional angles of the backbone obtained from

FIGURE 7. Torsional angles distribution obtainedfrom MC/1 (a) and MC/2 (b) simulations for the sevenresidues contained in the repeat of one helix of PMP.The population analysis has been performed consideringthe two torsional angles of each side chain as a pair.The three categories considered for each residue, in theorder displayed in the figure from left to right, are:trans-trans, trans-gauche+, and the other sequences.

JOURNAL OF COMPUTATIONAL CHEMISTRY 167

LEÓN ET AL.

TABLE II.Torsional Angles of the Backbone (in Degrees) andFiber Period of the Helix (in Å) for the 7/2 Helix ofPMP Obtained from the CLS Refinementsa and theMC/2 Simulation.

Number CLS MC/2

—C(HR)—C(H2)—C(HR)—CH2— −72 −77—C(H2)—C(HR)—CH2—C(HR)— 165 180Fiber period (c) 13.68 14.24

a From ref. 31.

CLS refinements with those obtained from theMC/2 simulation. As can be seen, the torsionalangles —C(HR)—C(H2)—C(HR)—CH2— andC(H2)—C(HR)—CH2—C(HR) resulting fromthe latter method differ by about 5◦ and 15◦,respectively, from those determined by the formerone. These differences involve a change in thefiber period length, i.e., the distance between thenearest equivalent residues of the 7/2 helix, of0.56 Å (Table II). This small difference, i.e., about4%, is not in contradiction with the experimentalX-ray diffraction diagram and with the measureddensity. Accordingly, it can be concluded thatthe stabilization in the packing of the PMP arisesfrom small changes in the torsional angles of thebackbone. It is worth noting that MC simulationsvarying the fiber period length but fixing the helixsymmetry are extremely useful to study systemswith poor X-ray diffraction data.

NPT Simulations of the Comb-LikePoly(α-octyl-β-L-aspartate)

Poly(α-alkyl-β-L-aspartate)s are stereoregularnylon 3 derivatives bearing an alkoxycarbonylgroup attached to the backbone β-carbon atom ofthe repeating unit.33 – 35

These polymers are reported to adopt helical con-formations similar to the α-helix typical of poly(α-amino acid)s. The most frequent crystal structureobserved in this family of compounds consists of a

monoclinic lattice composed of right-handed 13/4-helices arranged to be antiparallel. Such helices arestabilized by intramolecular hydrogen bonds set be-tween the i and i+ 3 amide groups.

The structure of poly(α-octyl-β-L-aspartate), ab-breviated PAALA-8, in the solid state was re-cently examined by IR dichroism, NMR, and X-raydiffraction.16, 17 The observed crystal structure wasinterpreted as an arrangement of sheets made of13/4 helices with side chains occupying the inter-sheet space (Fig. 1). The basic unit cell parametersobserved for this layered sructure are a = 12.3 Å,b = 18.0 Å, and c = 19.9 Å. It should be notedthat the distance between successive sheets, i.e.,18 Å, is not consistent with a fully extended con-formation for the side chain. Accordingly, the sidechains must be folded, especially those filling theinterchain space in the layers. On the other hand,a simple model of this biphasic structure was ob-tained by energy minimization.17

NPT simulations of PAALA-8 were performedwith the MCDP program. The model studied con-sisted of four independent helices of PAALA-8,each of them containing 13 residues. Helices werearranged antiparallel with respect to each other,in accordance with experimental data obtained forother poly(β-L-aspartate)s.34, 35 The helices werepacked in an orthogonal simulation box with cellparameters ao = 28.0 Å, bo = 39.0 Å, and co = 19.9 Å.It is worth noting that the starting a and b edgesare larger than those derived from the experimen-tal measures, i.e., a = 2 ∗ 12.3 Å and b = 2 ∗ 18.0.Thus, our aim is to achieve a correct description ofthe NPT behavior, which is a rather stringent test forany theoretical model.

The degrees of freedom in the simulations werethe torsional angles of the side chains and the lengthof the a and b edges of the cell. The torsional anglesof the 13/4 helix were kept fixed at the values re-ported in the literature,34, 35 which is consistent withthe experimental observation that this helical con-formation is even retained at high temperatures16

and in solution.36 Helix centers, defined as the inter-sections of helix axes with the xy-plane, were fixedrelative to the unit cell. In the CBMC algorithm, atotal of Ns = 6 torsional angles were used to samplethe torsional space for side chains. In addition to CBmoves, a small fraction (20%) of Metropolis moveswere also used. The simulations were performed atT = 298 K.

Figure 8a shows the energy as a function ofthe number of MC steps. As can be seen, the sys-tem was equilibrated rapidly Additionally, a blockanalysis indicates that the structures generated by

168 VOL. 22, NO. 2

MCDP

FIGURE 8. Energy (a) and cell parameters (b) ofPAALA-8 as a function of the number of steps.

MCDP are almost completely decorrelated after asfew as 10 cycles of moves for PAALA-8. Pilot cal-culations were performed in other poly(α-alkyl-β-L-aspartates) with longer alkyl side chains, denotedPAALA-n, where n is number of carbon atoms in thealkyl side group. Block size for complete decorrela-tion increases to about 1000 cycles for the longestside chains tested (PAALA-22), still many orders ofmagnitude more efficient than Metropolis. For evenlonger chains, Concerted Rotation5 moves will berequired to efficiently sample the torsional degreesof freedom closest to the main chain helix.

Figure 8b presents the values of the cell parame-ters a and b obtained from MCDP simulations. It isworth noting that the cell parameter a converges tothe a value, i.e., about 13.3 Å, close to experimen-tal one. On the other hand, b evolves spontaneouslyfrom the initially too large value towards 18.1 Å,which is also in excellent agreement with the exper-imental value. According to these results, it can bestated that MCDP simulations are able to reproducethe NPT behavior of PAALA-8.

Figure 9 shows a typical configuration forPAALA-8 obtained along the MCDP simulation,where the predominantly disordered state of the

FIGURE 9. Representative microstructure of thePAALA-8 obtained along the MCDP simulation.

interlayer region is qualitatively evident. To get adeeper insight in the side chain structure, the corre-sponding torsional angles for the alkyl side chainsof the 13 residues contained in the helix repeat ofPAALA-8 were analyzed. Figure 10 shows a popu-lation analysis for each torsional angle, the confor-mations being grouped in four categories: gauche+,

FIGURE 10. Torsional angles distribution for the 13residues contained in the helix repeat of PAALA-8. Thepopulation analysis of the torsional angle associated toeach bond in the side chain is specified. The fourcategories considered for each bond, in the orderdisplayed in the figure from left to right, are: gauche+,trans, gauche−, and the remaining conformers.

JOURNAL OF COMPUTATIONAL CHEMISTRY 169

LEÓN ET AL.

trans, gauche−, and the remaining conformers, i.e.,skew+ and skew−.

As it can be seen, in general, the trans is the pre-ferred conformation for the alkyl side chains. Thus,for torsions associated with the last 6 C—C bondsthis conformation is clearly the most populatedone. However, in all these bonds the populationof folded conformations is not negligible, whichprecludes the crystallization of alkyl chains in theinterhelical region. These results are fully consistentwith the partial order along the y-axis displayed inFigure 9 and experimentally observed for PAALA-8.

On the other hand, Figure 10 reveals that the tor-sional angles associated to the ester group and thebond number 2 of the alkyl chain prefer the foldedconformations. This is not a surprising result, be-cause these torsional angles are involved in the sidechain bending region, where the presence of foldedconformations is required to orient the side chainperpendicular to the planes of the helix sheets.

According to these results it can be concludedthat the MCDP strategy provides an excellent de-scription of both the NPT behavior and the confor-mational properties of PAALA-8.

Discussion and Conclusions

MCDP is a program that permits the simula-tion of dense polymers with ordered structures. Thestrategy combines two algorithms previously de-veloped by other authors: the CBMC algorithm1, 3, 4

and PCSP method.19 The CBMC method was usedby different authors to study a number of chemicalproblems in polymeric systems.3, 4, 7, 20, 37 – 41 How-ever, the CBMC algorithm was never implementedfor studying polymers with ordered structure, likefor instance comb-like polymers. The PCSP isa method explicitly developed to investigate thepacking of crystalline polymers. Thus, this methodallows one to do all the transformations requiredto study crystal structures. However, it is efficientonly for polymers with rigid conformations, likefor instance conventional nylons.22 – 24 The MCDPmethod is the result of combining both CBMC andPCSP.

MCDP is a useful tool for studying both the pack-ing of complex crystalline polymers and the struc-ture of comb-like polymers. Thus, MCDP includesthe full functionality of PCSP, because the conforma-tional flexibility of both the backbone and the sidechain is included. On the other hand, MCDP allowsatomistic simulations of comb-like polymers that ei-ther become impractical or lead to poor results with

other methods like molecular dynamics or energyminimization techniques.17

A problem that has to be considered is the size ofthe systems that can be simulated with MCDP. Toobtain the maximum computational efficiency, theprogram has been optimized and parallelized. Thespeedup reached on shared-memory multiproces-sors indicates that the program can be used for verylarge systems with a low computational cost.

The MCDP program can be used to investigatea wide variety of chemical properties in polymersystems related with their structure. For instance,the MCDP program has been recently used tostudy both the occupied space and the solubilityof small penetrants in comb-like poly(α-alkyl-β-L-aspartate)s with 12, 14, and 16 carbon atoms inthe alkyl group. Thus, the required structural in-formation, which is provided by MCDP, has beenextracted using an interface program, which allowsone to measure the volume not occupied by theatoms of the system and the excess of chemical po-tential for a given penetrant.18

Acknowledgments

M.L. gratefully acknowledges partial financialsupport from an Iberdola grant. S.L. acknowledgesthe support of the Ministry of Education of Spainfor the award of a scholarship. The authors thankto CEPBA (CEntre de Parallelisme de BArcelona)for computational facilities. We wish to thank thesupporting group of CEPBA for their help in theparallelization of the program.

References

1. Leontidis, E.; Forrest, B.; Widmann, A. H.; Suter, U. W.J Chem Soc Faraday Trans 1995, 91, 2355.

2. Binder, K. In Computational Modeling of Polymers; Bicer-ano, J., Ed.; Marcel Dekker: New York, 1992.

3. de Pablo, J. J.; Laso, M.; Suter, U. W. J Chem Phys 1992, 96,2395.

4. Siepmann, J. I.; Frenkel, D. Mol Phys 1992, 75, 59.5. Dodd, L. R.; Boone, T. D.; Theodorou, D. N. Mol Phys 1993,

78, 961.6. Krishna Pant, P. V.; Han, J.; Smith, G. D.; Boyd, R. H. J Chem

Phys 1993, 99, 597.7. Leontidis, E.; de Pablo, J. J.; Laso, M.; Suter, U. W. Adv

Polym Sci 1994, 116, 283.8. Olaj, O. F.; Lantschbauer, W. Macromol Chem Rapid Com-

mun 1987, 3, 847.

9. Carmesin, I.; Kremer, K. Macromolecules 1992, 21, 2189.

170 VOL. 22, NO. 2

MCDP

10. Zifferer, G.; Neubauer, B.; Olaj, O. F. J Chem Phys 2000, 112,428.

11. Watanabe, J.; Goto, M.; Nagase, T. Macromolecules 1987, 20,298.

12. Watanabe, J.; Ono, H.; Uematsu, I.; Abe, A. Macromolecules1985, 18, 2141.

13. Watanabe, J.; Takashina, Y. Polym J 1992, 24, 709.

14. Vierheller, T. R.; Foster, M. D.; Schmidt, A.; Mathauer, K.;Knoll, W.; Wegner, G.; Satija, S.; Majkrzak, C. F. Macromole-cules 1994, 27, 6893.

15. Poche, D. S.; Daly, W. H.; Russo, P. S. Macromolecules 1995,28, 6745.

16. López–Carrasquero, F.; Montserrat, S.; Martínez de Ilar-duya, A.; Muñoz–Guerra, S. Macromolecules 1995, 28, 5535.

17. Navas, J. J.; Alemán, C.; López–Carrasquero, F.; Muñoz–Guerra, S. Polymer 1997, 38, 3477.

18. Zanuy, D.; Alemán, C.; López–Carrasquero, F.; Báez, M. E.;García–Alvarez, M.; Laso, M.; Muñoz–Guerra, S., submittedwork.

19. León, S.; Navas, J. J.; Alemán, C. Polymer 1999, 40, 7351.

20. Siepmann, J. I.; McDonald, I. R. Mol Phys 1992, 75, 255.

21. Martin, M. M.; Siepmann, J. I. J Phys Chem B 1999, 103, 4508.

22. Bermúdez, M.; León, S.; Alemán, C.; Muñoz–Guerra, S.Macromol Chem Phys 1999, 200, 2065.

23. Bermúdez, M.; León, S.; Alemán, C.; Muñoz–Guerra, S. JPolym Sci Polym Phys 2000, 38, 41.

24. Bermúdez, M.; León, S.; Alemán, C.; Muñoz–Guerra, S.Polymer 2000, 41, 8961.

25. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller,A. H.; Teller, E. J Chem Phys 1953, 21, 1087.

26. Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids;Clarendon Press: Oxford, 1987.

27. Navas, J. J.; Alemán, C.; Muñoz–Guerra, S. Polymer 1996,37, 2589.

28. Voight–Martin, I. G.; Simon, P.; Yan, D.; Yakimansky, A.;Bauer, S.; Ringsdart, H. Macromolecules 1995, 28, 243.

29. Brisse, F. J Electron Microsc Technol 1989, 11, 272.30. Ferro, D. R.; Brückner, S. Macromolecules 1989, 22, 2359.31. Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A.

J Comput Chem 1986, 7, 230.32. Kusanagi, H.; Takase, M.; Chatani, Y.; Tadokoro, H. J Polym

Sci Polym Phys Ed 1978, 16, 131.33. Fernández–Santín, J. M.; Aymamí, J.; Rodríguez–Galán, A.;

Muñoz–Guerra, S.; Subirana, J. A. Nature 1984, 311, 53.34. Navas, J. J.; Alemán, C.; López–Carrasquero, F.; Muñoz–

Guerra, S. Macromolecules 1995, 28, 4487.35. García–Alvarez, M.; Martínez de Ilarduya, A.; León, S.;

Alemán, C.; Muñoz–Guerra, S. J Phys Chem A 1997, 101,4215.

36. Martínez de Ilarduya, A.; García–Alvarez, M.; Alemán, C.;López–Carrasquero, F.; Muñoz–Guerra, S. Macromolecules1999, 32, 3257.

37. Widmann, A. H.; Laso, M.; Suter, U. W. J Chem Phys 1995,102, 5761.

38. de Pablo, J. J.; Laso, M.; Siepmann, J. I.; Suter, U. W. Mol Phys1993, 80, 55.

39. de Pablo, J. J.; Laso, M.; Suter, U. W. J Chem Phys 1992, 96,6157.

40. de Pablo, J. J.; Laso, M.; Suter, U. W. Macromolecules 1993,26, 6180.

41. Siepmann, J. I.; McDonald, I. R. Langmuir 1993, 9, 2351.

JOURNAL OF COMPUTATIONAL CHEMISTRY 171