Accelerating electrostatic pair methods on graphical processing units to study molecules in supercritical carbon dioxide †

Traditionally, electrostatic interactions are modelled using Ewald techniques, which provide a good approximation, but are poorly suited to GPU architectures. We use the GPU versions of the LAMMPS MD package to implement and assess the Wolf summation method. We compute transport and structural properties of pure carbon dioxide and mixtures of carbon dioxide with either methane or di ﬂ uoromethane. The di ﬀ usion of pure carbon dioxide is indistinguishable when using the Wolf summation method instead of PPPM on GPUs. The optimum value of the potential damping parameter, a , is 0.075. We observe a decrease in accuracy when the system polarity increases, yet the method is robust for mildly polar systems. We anticipate the method can be used for a number of techniques, and applied to a variety of systems. Substitution of PPPM can yield a two-fold decrease in the wall-clock time.


Introduction
2][3][4] Maximizing computational throughput on graphical processing units (GPUs)/hybrid architectures is of great interest.Decreasing the wall-clock time for each MD time step allows longer timescales to be sampled, giving greater condence in the convergence of ensemble averages.GPUs have been widely adopted in various computational disciplines, due to the highly parallel design, low power consumption (FLOPS/Watt) and commodity cost (£/FLOPs).The GPU architecture is optimised for processing massive amounts of parallel calculations, and so has many more cores than a CPU, but at the sacrice of memory capabilities.
Electrostatic interactions can be divided into rst-order effects, which comprise point charge interactions decaying reciprocally with respect to intermolecular distance and higher-order interactions that decay more rapidly.
Electrostatic interactions can be modelled explicitly by various methods, including: cut-off truncation, [5][6][7] switched/shied cut-off truncation, 8 Ewald 9 and its mesh derivatives, particle mesh Ewald (PME), 10 particle-particle particle mesh (PPPM) 11 and smoothed particle mesh Ewald (SPME). 12The Ewald summation provides a good approximation for the electrostatic energy, as the algorithm accounts for periodicity of the simulation domain, but it scales as O(N 3/2 ).An approximation that is frequently adopted is PPPM, due to its superior scaling, O(NlogN), over the Ewald summation.SPME, 13 PPPM 14 and multi-level summation 15 have been implemented on the GPU, yet the requirement for fast Fourier transforms (FFTs) intrinsic to Ewald methods, like SPME and PPPM, reduces the parallelism 14 and leads to poor scaling on GPUs.Pair-wise algorithms derived from cut-off techniques show superior scaling O(N) and greater effectiveness on GPUs.
Cut-off techniques can reproduce experimental quantities, for instance, the Madelung constant.Wolf et al. 16 reported that the electrostatic (Coulombic) potential for condensed systems was effectively short-ranged, and the energies are in agreement with Ewald methods when the cut-off sphere reached neutrality.This was observed whilst calculating the Madelung constant for rocksalt (NaCl), where the relationship between truncation distance and electrostatic energy was investigated.When the cut-off sphere was truncated at a distance to achieve charge neutrality, the electrostatic energy was signicantly closer to the Madelung energy.The repeating lattice structure of NaCl is well-suited for treatment with the Wolf method, but the supercritical phase of polar solutes can be more difficult to model as the atomic partial charges within the cut-off sphere are dynamic.The resulting shied Coulomb potential achieves charge neutrality by projecting each atom charge, q j , from q i onto the edge of the solvation sphere.Therefore, every jth atom in the solvation sphere of atom i has a charge of equal but opposite sign set at the cut-off (R c ).This results in an articially neutral solvation sphere for every ith atom, which effectively makes the system charge neutral.
where q i and q j are point charges, r ij is the intermolecular distance between atoms i and j, a is the damping coefficient in ÅÀ1 and R c is the cut-off distance.To obtain accurate electrostatic contributions, a damping function is applied.The electrostatic energy decay oscillates around a rate of 1/R c .Introducing the damping function quickly attens the oscillations as the cut-off increases, effectively determining how fast the complementary error function falls from unity at r ij ¼ 0, to zero at the cut-off. 16The damping function adopted is the complementary error function, as it has a close connection to the Ewald sum. 17The coefficient of the error function, a, denotes the rate at which convergence is achieved.A large a value will converge the energy rapidly using a short cut-off, but the errors can be larger.A smaller a leads to less contamination of the potential, but the sum will uctuate more rapidly.Assigning a value of 0 to a, results in the truncated shiedforce (SF) potential.This SF calculation is faster than the Wolf method, but the selection of a can enhance the accuracy of electrostatic forces and energies by optimizing agreement to Ewald methods. 18imple cut-off based methods are unreliable for computing the forces, as the potential truncates abruptly at the cut-off, which causes forces to be undened if the molecule lies at the cut-off boundary.The reliance upon Ewald methods instead of efficient cut-off methods has been discussed by Fennel and Gezelter, 17 who compared the accuracy of damped shied-force (DSF) and shied potential cut-off methods for water, argon in water, and NaCl.Benchmarking and validation of the DSF potential was applied to polyelectrolyte brushes on GPUs; speedups were achieved of a factor between 1.1 and 3.9, depending on the system and size of the cut-off with respect to PPPM. 19Group based cut-offs, where the cutoff is dynamically allocated to ensure entire molecules are within the solvation shell, were investigated, 17 but the energies deviated signicantly from those obtained using PME.
In this work we assess the applicability of the Wolf method to model electrostatic interactions for various systems in the supercritical region of carbon dioxide.Supercritical carbon dioxide (scCO 2 ) is an attractive 'green solvent' used in many industrial processes, such as caffeine extraction, 20 polymer solvation 21 and synthesis, 22,23 enzyme catalysis 24,25 and stabilization 26,27 and transition metal catalysis. 28,29Although carbon dioxide does not possess a dipole, it has a signicant quadrupole 30 (13.4 Â 10 À40 Cm 2 ), which means that electrostatic interactions are an important component of interactions involving scCO 2 .Su and Maroncelli 31 observed that neglecting electrostatic interactions, i.e., considering only Lennard-Jones interactions, led to a systematic 14% error in the solvation free energies of polymer-scCO 2 systems.This observation is attributed to quadrupole-dipole and quadrupole-quadrupole interactions, 31 which are inherently modelled in all point charge models.A nonpolar uid should be well suited for treatment using the Wolf method, as quadrupolar interactions decay more rapidly than dipolar ones.This bodes well for the efficient electrostatic treatment of systems solvated by carbon dioxide by utilizing the Wolf method on GPU architectures.We also investigate the computational cost of the DSF method, with respect to the Wolf method.
Fluorous polymers are well known to dissolve in scCO 2 , 32,33 which has been partially explained by 19 F-NMR experiments, which suggested a number of specic interactions between carbon dioxide and the uorous solute that increase the solubility. 34Many biomolecular systems are stable in the presence of scCO 2 , but this is dependent on the species, the water content, and experimental conditions. 24Protein stability can be observed when scCO 2 solvates the hydrophobic residues, and water solvates the polar/hydrophilic residues. 26Water possesses a strong dipole moment, which will make it less amenable to treatment with the Wolf method.We aim to follow this work with an investigation for protein systems solvated by water and carbon dioxide.
To assess the applicability of the Wolf method to systems containing carbon dioxide, important physical quantities, such as PVT relationships and diffusion coefficients, can be calculated and compared with the values calculated using PPPM.The applicability of the Wolf method to study methane plus carbon dioxide gas hydrates has been reported, and the results show good agreement between lattice sum and reaction eld methods. 35Analysis of the convergence behaviour of the Wolf method by Angoshtari and Yavari 36 shows the method to be robust, but convergence can be problematic if poor choices are made for the cutoff or a values.In our study, we investigate the effect of increasing the polarity of the uid by incorporating diuoromethane molecules, which possess a strong dipole.For low polarity systems the Wolf method should be well suited, but as the polarity increases the effect of long-range dipole interactions will become important and we assess the point where Ewald techniques will become necessary.

Method
LAMMPS is a multi-purpose MD code, 37 which is widely used in the elds of atomistic, coarse-grained and mesoscopic simulations.The code can process massive numbers of particles per simulation, which it achieves using an optimized spatial decomposition technique.The GPU optimized version, written in CUDA C, can run over multiple GPUs, either in conjunction with the CPU or entirely on the GPU.We have implemented the Wolf method into LAMMPS-CUDA the GPU-exclusive version of LAMMPS and LAMMPS-GPU the CPU/GPU implementation, as a new potential incorporating Lennard-Jones and electrostatic interactions.LAMMPS-CUDA was written exclusively for use with CUDA, but LAMMPS-GPU can be compiled and run on AMD GPUs and other applicable accelerators using OpenCL.The new pair style is implemented as a double potential function in LAMMPS as lj/charmm/coul/wolf for CPU, lj/charm/coul/ wolf/cuda for LAMMPS-CUDA and lj/charm/coul/wolf/gpu for LAMMPS-GPU.The LAMMPS-CUDA version is solely GPU based, with the exception of le I/O and pre/ post simulation setup.The GPU version of LAMMPS utilizes the GPU for the force and/or neighbour list generation, whilst all other operations and I/O are performed on the CPU.The GPU version allows for n CPUs to be used per GPU, whilst the LAMMPS-CUDA version only allows for one CPU per GPU.An important difference between packages is that the GPU version on LAMMPS uses the CPU to calculate FFTs for PPPM, whilst the LAMMPS-CUDA version performs the calculation on the GPU.We have also implemented the DSF implementations, lj/ charm/coul/dsf/cuda and lj/charm/coul/dsf/gpu, for benchmarking purposes.We used the CHARMM 38 Lennard-Jones potential in the AB form with a switching function in conjunction with the Wolf method.A cut-off of 2.5s was used for the Lennard-Jones potential, and tapered to zero at 2.65s, where s was calculated Table 1 EPM2 force field parameters for carbon dioxide 39 and Palmer and Anchell parameters for methane and difluoromethane. 40(A 1/12 ) and (B 1/6 ) have units ((kcal mol) 1/12 Å) and ((kcal mol) This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.
from A 1/12 in Table 1.We used the rigid EPM2 atomistic force eld 39 to represent carbon dioxide and parameters derived by Palmer and Anchell 40 to represent Lennard-Jones and charge interactions for methane and diuoromethane.We utilized the LAMMPS soware to simulate scCO 2 at several densities, comparing the results obtained by using PPPM or the Wolf method.We selected a tolerance setting of 0.0001 for the calculations involving PPPM that enables the root mean square error to be within a factor of 10 000 of the reference force, which is calculated analytically for short-range interactions 41 and in k-space. 42The optimum value of a was investigated for binary mixtures of carbon dioxide and methane or diuoromethane.The criterion for selecting a is the level of agreement between the Coulombic energies computed by PPPM and Wolf.To consider the effects of neglecting only the periodic long-range Coulombic effects, we treat the cut-off as being half the length of the simulation box.We also investigate the level of approximation that arises from use of a quarter box cut-off in conjunction with the Wolf method.

Benchmarking GPU electrostatics
Two boxes of 10 000 and 50 000 carbon dioxide molecules were constructed with a molar density of 0.01 mol cm À3 and minimized for 1000 iterations using the conjugate gradient method, to an energy tolerance of 1 Â 10 À6 kcal mol À1 .Both systems were heated from 0 K to 308.2 K using the Nosè-Hoover thermostat 43 with a damping time of 500 fs for 2 ns, followed by 5 ns of equilibration at 308.2 K.We evaluated the performance of the Wolf method on two GPU architectures, Tesla (Tesla C1060) and Kepler (Tesla K10).The K10, released in 2012, is designed for high throughput calculations performed in single precision.The number of giga oatingpoint operations per second (GFLOPs) is 4577 in single precision, compared to 933 for the C1060.This improvement was achieved in part by increasing the number of cores in the streaming multiprocessor from 8 (Tesla), 32 (Fermi) to 192 (Kepler).Memory bandwidth on the GPU to global memory is 103 GB s À1 for the C1060, but 320 GB s À1 for Kepler K10.We include an eight core Intel Xeon (E5-2609) CPU for benchmarking purposes that has an approximate 77 GFLOPs with 34 GB s À1 bandwidth to RAM.Calculations were performed using single-precision on two separate nodes, both comprise an eight core Xeon CPU and either a Kepler K10 or two Tesla C1060.Multiple GPU tests were only carried out using the C1060, which yielded an almost two-fold increase in performance, which corresponds well to the linear dependency.We use the FFTW 3.3.1 library in single-precision for all k-space calculations of PPPM.We considered two modes for calculations using LAMMPS-GPU, one where all force and neighbour calculations are performed on the GPU, and the other where the force and neighbour calculations are dynamically assigned between CPU and GPU.We used the CUDA 5.0 Tokito (GPU driver v. 304.54), which has produced a noticeable improvement in performance over CUDA 4.0.A new feature, CUDA dynamic parallelism, allows kernels running on the GPU to spawn more grids and to continue to generate work depending on the calculation. 44This feature has not been incorporated in our study.

Pure carbon dioxide
The 10 000 carbon dioxide molecule system was used to generate densities (box lengths) of 0.001 (255.1 Å), 0.002 (202.5 Å), 0.004 (160.72 Å), 0.005 (149.2Å), 0.01 (118.42Å) and 0.02 (94.0 Å) mol cm À3 , all within the supercritical region at 308.2 K.Each box was minimized for 1000 iterations using the conjugate gradient method, with an energy tolerance of 1 Â 10 À6 kcal mol À1 .The system was heated in the NVT ensemble using the Nosè-Hoover thermostat 43 with cubic periodic boundary conditions from 0 K to 308.2 K for 1 ns using a 1 fs time step.This was followed by 1 ns of equilibration.Integration was performed using the timereversible velocity Verlet algorithm. 45The procedure was repeated nine times to generate different equilibrated congurations for the purposes of error analysis, and the resulting standard deviation is noted in the error bars.The mean square displacement (MSD) of the centre of mass of carbon dioxide was obtained by calculating the gradient of the linear portion of the relationship between MSD and lag time over 10 ns production dynamics in the NVT ensemble.The pressures were obtained every 50 fs, and averaged over 10 ns to calculate the PVT relationship.

Carbon dioxide and solute
To study the effects of polar and non-polar solutes, we considered different quantities of methane and diuoromethane molecules.Five boxes of 10 000 carbon dioxide molecules were constructed with mole fractions c (solute) ¼ 0.0001, 0.001, 0.01, 0.09 and 0.5, where the solute was either methane or diuoromethane.These systems correspond to box lengths ( Å) of 95.01, 95.99, 97.49, 99.27 and 125.31 for CH 4 /CO 2 and 81.64, 81.78, 81.83, 84.44 and 105.16 for CH 2 F 2 /CO 2 .Each system was minimized and heated to 308.2 K using the same procedure as above.The pressure was equilibrated to 80 atmospheres (corresponding to a molar density of $0.02 mol cm À3 ), using the Berendsen barostat 46 with a 1 ps damping time for 20 ns with a 1 fs time step.The procedure was repeated nine times to generate different equilibrated congurations for the purposes of error analysis, and the resulting standard deviation is noted in the error bars.The MSDs of carbon dioxide and solutes were obtained from the linear portion of the relationship over a production run of 10 ns.We compare the errors in Coulombic energy with respect to a for c (solute) ¼ 0.09 with PPPM and Wolf, which we performed by decomposing the energy into group contributions.We investigate the relationship between system polarity and diffusion coefficients, and the total Coulombic energy.
One method of quantifying the transport properties of a system is the MSD, and thus the diffusion coefficient using the Einstein relationship.
where D is the macroscopic diffusion coefficient, r i (t) is the position of the centre of mass at time t, r i (0) is the initial position of the centre of mass, t lag is the lag time, and h[r i (t) À r i (0)] 2 i is the ensemble averaged MSD.The diffusion coefficient is calculated from the slope of the MSD against lag time, 47 which is a measure of the atomic displacements through time.
The pressure of a system can be calculated using the virial equation (below), and the PVT relationship has an inuence on the diffusion coefficients. 48The long-range part of PPPM has a different contribution to the pressure. 49here P is the pressure of the system, V is the volume, N is the number of particles, x is the dimensionality, K B is the Boltzmann constant and the term in brackets is the total intermolecular force multiplied by the interaction distances.

Benchmarking
To compare the efficiency of the Wolf method, we measured the computational throughput of the electrostatic routines using the 10 000 and 50 000 molecule systems.The timings of the implementations are calculated with respect to the number of neighbours in the cut-off sphere.Fig. 1(a) and 1(b) compare the throughput for different architectures, and show the acceleration gained when substituting PPPM for Wolf using LAMMPS-CUDA.The computational throughput is almost twice as high with the Wolf treatment than with PPPM.The Wolf method is slower than PPPM on the CPU, which is unexpected, but the CPU is well suited to handle FFTs that feature in PPPM.This is due to faster memory accesses between cache and compute cores, and the less parallel nature of the algorithm.The DSF method is marginally slower than PPPM in LAMMPS-CUDA, but is marginally faster than the Wolf method for the LAMMPS-GPU implementation.The difference between pair styles using LAMMPS-GPU is much less pronounced than LAMMPS-CUDA, where the acceleration gained by using the Wolf or DSF method is approximately 1.5 times faster.The difference between dynamic partitioning and exclusive GPU computation is small for the cut-off methods, and the dynamic mode is marginally quicker.For PPPM the opposite is observed, which could be due to the CPU already being used for the k-space calculations.The Kepler K10 GPU shows the greatest acceleration, which bodes well for future GPU releases.

Optimum alpha and cut-off parameters
We consider rst the accuracy of the electrostatic energy and its convergence as a function of a and cut-off distance.We then turn to the properties computed from the simulation, considering structural, dynamic and thermodynamic aspects.The simulations have all been executed in the supercritical region of carbon dioxide, and therefore each molecule has a large number of neighbours in the solvation sphere.We compare the effects of selecting a values between zero and 0.3 in intervals of 0.025.We decompose the energy interactions into pairwise group contributions to the total Coulombic energy for CH 4 in carbon dioxide, and CH 2 F 2 in carbon dioxide.The relationship between a and the accuracy of the Wolf summation method is shown in Fig. 2. We observe for half box cut-offs the best agreement for Wolf is at low values of a, where for both the polar and non-polar systems the optimum value for a is This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.
0.075.For non-polar systems increasing a beyond 0.2 results in an $2-3% error (with respect to the PPPM Coulombic energy), whilst for the polar system the average error is $20%.The non-polar CH 4 /CO 2 system resulted in the lowest average errors of 0.05% for CO 2 -CO 2 interactions, 0.39% for CO 2 -CH 4 and 0.34% for CH 4 -CH 4 .We observe a similar trend for the polar system of CO 2 /CH 2 F 2 with increased average errors of 0.44% for CO 2 -CO 2 , 0.47% for CO 2 -CH 2 F 2 and 0.67% for CH 2 F 2 -CH 2 F 2 .At a ¼ 0.075 the greatest variance in the non-polar system is s 2 ¼ 0.2 for CH 4 -CH 4 , and for the polar system the maximum variance is s A quarter box cut-off gives greater errors in the Coulombic energies, but follows the same trends as the half box cut-off.For the non-polar system the optimum value for a is 0.05, which is seen for all pair interactions of methane and carbon dioxide.
The average error increases to 0.5% for CO 2 -CO 2 and to $0.75% for energies involving interactions with methane.The errors vary more when used with the polar system over a wide range of a, but at a ¼ 0.075 the error for all interactions is below 2%.The error increases sharply for the polar system with varying a.For interactions between CH 2 F 2 -CH 2 F 2 the variance is higher than for methane, which can be attributed to the high polarity of both substituents.

Pure carbon dioxide
Diffusion coefficients of scCO 2 calculated from MD simulations using the Ewald sum 50 have been previously compared with experimental results. 51Our simulations show that the Wolf method gives diffusion coefficients comparable to that of PPPM (Fig. 3b).Calculations were concurrently run on the CPU using the Wolf method, and the energies were within three decimal places.Simulations were also performed at 323.2 K, which show good agreement between PPPM and the Wolf method.† Low density boxes using the Wolf method show the best agreement, but all the densities investigated are within the bounds of error of PPPM.Pressures obtained over a range of densities coincide well with an equation of state, 52 for both the Wolf and PPPM implementations (Fig. 3a).Both electrostatic methods capture the PVT relationship properties at high and low densities, for both temperatures studied.

Binary mixture of carbon dioxide with diuoromethane or methane
The Coulombic energy for the polar systems is compared for different mole fractions of diuoromethane using the Wolf method and PPPM.As the system becomes more polar, the error of the Wolf method with respect to the PPPM value becomes greater.Fig. 4 shows the total PPPM Coulombic energy, and results for a ¼ 0.2 and a ¼ 0.075.For the system containing 100 diuoromethane molecules, the error is within 0.2% and 1.8% for a ¼ 0.075 and a ¼ 0.2 respectively.As the composition tends towards a 1 : 1 ratio, the errors increase, which indicates, as anticipated, that the Wolf approach is not suitable when the system becomes highly polar.The errors for a ¼ 0.075 are reasonable (0.9%) but for a ¼ 0.2 the errors are 3%.We observe better agreement for methane.For the 1 : 1 binary mixture, the error in a ¼ 0.2 is 1.8% and a ¼ 0.075 is 0.15%.The total Coulombic energy for diuoromethane is more negative than for methane, which can be explained by NPT dynamics resulting in a lower volume box and therefore closer contacts.
Experimental results from O'Hern and Martin 51 indicate pure carbon dioxide has a diffusion coefficient of 5 (Â 10 8 m 2 s À1 ) at 308.2 K, 80 atmospheres.The nonpolar binary mixture of CH 4 /CO 2 has a diffusion coefficient of $5 (Â 10 8 m 2 s À1 ) for the inclusion of one methane molecule in 10 000 molecules of carbon dioxide (Fig. 5).As the number of methane molecules increases, the overall diffusion coefficients remain relatively constant.Diffusion coefficients for CH 4 and carbon dioxide are similar, with an average increase of $0.6 (Â 10 8 m 2 s À1 ) between carbon dioxide and CH 4 .As c (solute) increases, the accuracy of the solute diffusion coefficients increase, due to better averaging from more solute molecules.With the exception of c (solute) ¼ 0.0001 where averaging is poor, the agreement between Wolf and PPPM is good.We observe a reduction in diffusion coefficients for diuoromethane as the polarity increases.The values decay from 4.5 (Â 10 8 m 2 s À1 ), to 1.3 (Â 10 8 m 2 s À1 ), which can be attributed to favourable interactions between solute and solvent.Agreement between Wolf and PPPM is good, although there is an observable increase in the error in the diffusion coefficients for diuoromethane when the fraction of solute reaches 10%.This indicates the diuoromethane interacts strongly with carbon dioxide, thus limiting diffusion.
To characterise the interactions between diuoromethane and carbon dioxide, we calculate the radial distribution (RDF) between the centre of mass for CH 4 /CO 2 and CH 2 F 2 /CO 2 (Fig. 6) and the associated residence times and coordination numbers are shown in Table 2.
The RDF shows a larger density of carbon dioxide molecules in the diuoromethane mixture in the rst and second solvation shell compared to methane.The same trend was noticed in the RDF by Do et al., 53 where the rst solvation shell has a higher density for diuoromethane than for methane.The number of carbon dioxide molecules present in the rst and second solvation shells is $30% higher for diuoromethane, and carbon dioxide resides about four times longer compared to methane.This indicates that carbon dioxide has a higher affinity for diuoromethane.

Conclusions and discussion
The Wolf method shows good agreement with PPPM when modelling the electrostatic interactions of scCO 2 on GPUs for non-polar systems, whilst being approximately twice as fast.The choice of a is important, and may need to be investigated on a case by case basis to enable satisfactory agreement.For modelling carbon dioxide in the supercritical region it is advisable to use a halfbox cut-off and a low value for a.In this investigation, all values of a less than 0.15 produced errors less than 2% for non-polar interactions, whilst polar interactions require a to be less than 0.1.Upon increasing the polarity of the system, the  potential begins to degrade.Errors of the Wolf method with respect to PPPM are approximately 0.2% when considering a 100 : 1 mixture of carbon dioxide and diuoromethane with a ¼ 0.075.We observe a strong affinity of carbon dioxide to diuoromethane compared to methane, which can be seen by a decline in diffusion coefficients with increasing mole fraction of solute.Carbon dioxide resides about four times longer in the solvation sphere of diuoromethane compared to methane.We can conclude that the signicance of using the Wolf method on GPUs allows simulations to reach timescales twice as long as those run with PPPM, without signicant loss in accuracy for a carefully chosen value of a for non-polar and mildly polar systems.We aim to follow up the investigation with further analysis of solvent-solute interactions and the study of uorinated polymers, which have high solubilities in scCO 2 .We will be investigating the free energy changes of uorinated polymers, with an aim of further understanding the high affinity of uorous polymers for carbon dioxide.Many free energy methods require long timescales in order to reach convergence; utilizing the Wolf method for this purpose will help achieve this goal.

Fig. 1
Fig. 1 Computational throughput as a function of number of neighbours for (a,c) 10 000 CO 2 molecules (B-D) 50 000 CO 2 molecules.Fig.(a) and (B) were calculated using the LAMMPS-CUDA package, and (c) and (d) were calculated using the LAMMPS-GPU package.In (a) and (B), three electrostatic treatments were considered, Wolf (filled shapes), PPPM (unfilled shapes) and DSF (green circles) for different architectures.(Triangles ¼ Tesla C1060, squares ¼ Kepler K10, diamonds ¼ eight core Intel Xeon).In figures (c) and (d) the force/neighbour implementation is shown where the GPU is used exclusively (diamonds) and dynamically assigned between CPU/GPU (squares).PPPM is shown in black, the Wolf method is shown in blue, and the DSF force/neigh exclusively GPU is shown in green.

Fig. 2
Fig. 2 Percentage error in Coulombic energy between PPPM and Wolf with respect to a for CH 4 /CO 2 and CH 2 F 2 /CO 2 binary mixtures for c(solute) ¼ 0.1 at 308.2 K, 80 atmospheres.The results are decomposed into (a) (CH 2 F 2 ) CO 2 -CO 2 pair errors (B) (CH 2 F 2 ) CH 2 F 2 -CH 2 F 2 pair errors, (c) (CH 2 F 2 ) CO 2 -CH 2 F 2 pair errors, (d) (CH 4 ) CO 2 -CO 2 pair errors, (e) (CH 4 ) CH 4 -CH 4 pair errors and (f) (CH 4 ) CO 2 -CH 4 pair errors.The black squares and error bars indicate the average error and uncertainty for a half box cut-off, whilst the empty squares show the average error for a quarter box cut-off.

Fig. 3
Fig. 3 (a) Diffusion coefficients for pure carbon dioxide obtained for PPPM (empty red squares) and Wolf (filled blue squares), compared with experimental (black line) at 308.2 K. b) PVT relationship of pure carbon dioxide simulated using PPPM (empty red squares) or the Wolf method (filled blue squares) compared with experimental (black line) at 308.2 K.

Fig. 4
Fig. 4 Total Coulombic energy for methane and difluoromethane in carbon dioxide at 308.2 K, 80 atmospheres.The dashed red line indicates the total PPPM energy, whilst the green squares indicate a ¼ 0.2 and the blue squares indicate a ¼ 0.075.

Fig. 5
Fig. 5 Diffusion coefficients for (a) CH 4 in CO 2 /CH 4 (B) CO 2 in CO 2 /CH 4 (c) CH 2 F 2 in CO 2 / CH 2 F 2 and (d) CO 2 in CO 2 /CH 2 F 2 at 308.2 K, 80 atmospheres.Filled blue squares and error bars indicate diffusion coefficients obtained using PPPM, and empty red squares and error bars indicate diffusion coefficients obtained using the Wolf summation method.

Fig. 6
Fig. 6 Radial distribution function between the centre of mass of carbon dioxide and CH 4 (solid line), and CH 2 F 2 (dashed line) for c (solute) ¼ 0.1 at 308.2 K, 80 atmospheres.

Table 2
Residence times and coordination for the centre of mass of carbon dioxide in CH 4 and CH 2 F 2 for c (solute) ¼ 0.1 at 308.2 K, 80 atmospheres Faraday Discussions Paper 354 | Faraday Discuss., 2014, 169, 343-357 This journal is © The Royal Society of Chemistry 2014 Open Access Article.Published on 05 March 2014.Downloaded on 12/11/2018 3:08:52 PM.This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.