Results and Analysis

Using the the GCCBMC algorithm, several physical properites of polymeric systems were studied in this project. We studied the temperature and density effects on polymer conformation and the phase behavior of the polymer system. The efficiencies of the GCCBMC algorithm is also compared with an algorithm using reptation to move the polymers.
Our first concern is to make sure that the simulation is working properly and also to characterize the system. To do this, the radius of gyration of the polymer (Rg) is measured at various chainlengths (N_m) and temperature (T). Statistical theories of polymers⁴ predict that Rg ~ N_m^v , where v is a temperature - dependent coefficient due to competing intra-polymer and inter-polymer interactions. Theories, experiments, and other simulations have shown that at a theta temperature (T_theta), the polymer has an ideal random walk conformation and v = 0.5. For T >T_theta, the polymer has a self-avoiding walk conformation and v=0.6. For T < T_theta, the polymer has a collapsed conformation and v = 0.33. The system used to study the polymer conformations in this project has fixed volume, temperature, and only one polymer. The polymer is allowed configurational-bias moves. The results of the simulations done are shown in the following plot.

The exponent v is extracted from the simulations of a single polymer performing configuration moves in a NVT system of boxlength 100. By varying the number of monomer units (N_m) in the polymer from 4 to 64, the exponent is determined from curve fitting Rg vs N_m.(Click on the graph to get an enlarged view.)
Another physical trend that we can check against is that the monomer - monomer correlations between polymers should increase as the density of polymers increase. The inter-polymer density correlation is defined as

where V is the system volume, N is the number of particles in the system, and r_ij is the distance between particles i and j. The algorithm used to calculate g(r) from the simulation data can be found here. The inter-polymer density correlation function was calculated from the simulations done at T = 4, boxlength = 12. 1000 MC steps was collected after 2000 MC steps of equilibration.
From the graph, it can be concluded that the simulation does agree with the physical trend of increasing correlation with increasing density. Also, at a higher density of 0.22 (24 chains), a small peak at r = 1 can be observed. This indicates that polymers are starting to form structures, and this effect is more prounouced at higher densities, as shown in the following plot.(Click on the graph to get an enlarged view.)

This plot shows a NVT system of boxlength 8.6 at T=1.9 containing 8-monomer chains. The peaks in the density-correlation functions are now prounced and indicate liquid shells of polymers forming in the system.(Click on the graph to get an enlarged view.)

Our next concern is about the efficiency of the configurational - bias algorithm. Because polymers are molecules on chains, it is very difficult to move them around with random walk in a dense system using conventional Monte Carlo due to the low acceptance ratio. Algorithms such as reptation⁷, pivoting⁸, and configurational bias have been used by researchers to increase the acceptance ratio. In this project, we compared the efficiencies of the configurational bias method against that of the reptation method using Rg as the quantity of concern. Since efficiency ~ 1/ (error^2 * time), the error bars in Rg and computer time are determined from the simulation.
For this comparison, we chose a system with 50 polymers, each polymer consists of 8 monomers. The system is a cube with length 8.6 to a side, and the temperature is chosen as 1.9. The density of this system is 50*8/(8.6)^3 = 0.63. The system is allowed to equilibrate for 2000 Monte Carlo steps, then the data is sampled for 500 MC steps. The following plot shows that using configuartional-bias algothrim allows faster equilibration and also achieves better efficiency under these conditions.(Click on the graph to get an enlarged view.)

Comparison of Efficiencies: 1000 MC steps after 4000 steps equlibration
time Rg sigma(Rg) corr. time efficiency acc ratio (%)

Reptation 119.7 1.420 0.011 50 15.6 3.78

CBMC(k=8) 400.9 1.412 0.007 15 23.8 10.97

Comparison of Efficiencies: 1000 MC steps after 4000 steps equlibration
	time	Rg	sigma(Rg)	corr. time	efficiency	acc ratio (%)
Reptation	119.7	1.420	0.011	50	15.6	3.78
CBMC(k=8)	400.9	1.412	0.007	15	23.8	10.97

In a less dense system of 20 polymer chains, the following data shows that CBMC is much more efficient than reptation.(Click on the graph to get an enlarged view.)

Comparison of Efficiencies: 1000 MC steps after 2000 steps equilibration
time Rg sigma(Rg) corr. time efficiency acc. ratio(%)

Reptation 55.5 1.446 0.012 100 15.5 19.8

CBMC(k=8) 174.6 1.326 0.0045 30 42.4 33.9

Because we do not use a potential cutoff radius in the simulations, the CBMC method becomes less efficient as the number of particles increases. If a cutoff radius is employed, the CBMC method should be even more efficient relative to the reptation method as system density increases.

Having characterized the dependence of polymer conformation on temperature and determined that CBMC is more efficient than reptation. Our studies then focus on adsorption and phase coexistence curve of polymeric systems.

The adsorption curve of a polymer system is obtained by using a grand canonical ensemble and allow polymers from a reservoir to enter the simulation box. The initial configuration of the simulation box has boxlength 8.6, T=1.9, and contains 50 polymers of 8-monomer units. After 2000 MC steps for equilibration, 2000 more were run for data collection. The plot is shown as follows. The density in this plot is the density of polymer chains.(Click on the graph to get an enlarged view.)
Another important part of understanding the phase behavior of the polymer system is to find the phase coexistence curve if there is one. In this project, this was done by using the Gibbs ensemble. In the algorithm, this means making two simulation boxes of constant total volume, and allow the individual box volumes to change and the polymers inside these boxes to interchange. Initially, each box contains 50 polymers of 8-monomer units. After 6000 MC steps of equilibration, 6000 MC data are collected. The two systems separate into a higher density system (liquid) and a lower density system (vapor) as shown in the following plot.
From this plot, the critical temperature can be estimated as around T_c = 1.9 +/- 0.1, and the critical density is around_c = 0.08 +/- 0.002. (rho is defined in this plot as # of polymer chains / Volume.)(Click on the graph to get an enlarged view.)
Furthermore, we can plot the oder parameter_L -_G vs temperature T.(Click on the graph to get an enlarged view.)

We can extract the critical exponent b from _L -_G
~ (1-T/T_c)^b .If we fit the first 4 data points with T=1.9, we get b=0.34+/-0.025. If we fit the first 7 data points with T=1.8, we get b=0.44+/-0.08. We need more data points of temperature between 0.8 and 1.5 in order to extract a meaningful value.

Back to Project Main Page