research on lipid membranes 
 and associated proteins has broad importance 
 in understanding various biological processes such as lipid trafficking , 
 antibiotic resistance , drug transport , cell signaling , and lipid raft 
 formation .
molecular dynamics ( md ) simulation has become an important 
 technique to study such systems , as computational resources and empirical 
 force fields have advanced to allow for meaningful simulations of 
 lipid membranes to better understand their physical properties and 
 functions .
for example , md simulations have been used to probe details 
 of liquid - ordered and liquid - disordered phases and also the formation 
 of liquid - ordered rafts .
there has also been a recent focus 
 on developing accurate representations of cellular membranes based 
 on the known diverse lipid compositions in cells and organelles .
specifically , md simulations have been used to model the plasma 
 membrane at a coarse - grained level to fully atomistic 
 representations of inner and outer membranes of gram - negative bacteria .
the transport and ligand binding mechanisms for primary and secondary 
 active transporters and binding of peripheral membrane proteins to 
 the bilayer surface have also been studied using md simulation techniques .
important to md simulation is an accurate mathematical description 
 of the relationship of conformation to energy in and between molecules , 
 which is commonly referred to as an empirical force field ( ff ) . for 
 lipids , the all - atom additive charmm36 ( c36 ) ff was developed to accurately represent a wide range of experimental 
 bilayer properties using the constant particle number , pressure , and 
 temperature ( npt ) ensemble .
the c36 ff has been extended to a wide 
 range of lipids , e.g. , sphingolipids , phosphoinositides , cardiolipin , bacterial lipids , ester - modified 
 lipids , those with polyunsaturated chains , cholesterol , and 
 lipopolysaccharides .
recently , the amber community 
 has developed their own all - atom lipid ff ( lipid14 and slipids ) that is of comparable accuracy to that of the c36 
 ff .
one important caveat in these ffs and their usage is that they 
 are developed using approaches that are consistent within the family 
 of each ff , so that it is generally desirable to remain internally 
 consistent ( within the same ff family ) when studying heterogeneous 
 systems and avoid mixing the parameters of different ffs without careful 
 testing .
moreover , methods used for truncation of the nonbonded interaction 
 energy should be maintained as those used during optimization of the 
 ff to accurately represent these ffs during application studies .
the form of the potential energy function used in most biomacromolecule - based 
 ffs are highly similar and thus allows various ffs to be used in different 
 md simulation programs .
the c36 ff was developed by the charmm community 
 initially in the context of the charmm program with the treatment of the long - range nonbonded interactions 
 as follows .
the lennard - jones ( lj ) 6 - 12 ( i.e. , van der waals ) interactions 
 were smoothed over the range of 10 to 12  using the force - based 
 switching function , the particle mesh 
 ewald ( pme ) method was used for the long - range 
 electrostatic interactions , and 1 - 4 nonbonded interactions ( i.e. , 
 coulombic and lj interactions for atoms separated by three chemical 
 bonds ) are fully turned on ( i.e. , no scaling ) . the cutoff and
1 - 4 
 interaction scheme varies between the ffs and should be maintained 
 when different simulation programs are used .
this is particularly 
 important with lipids , which have been shown to be particularly sensitive 
 to the treatment of the truncation of the lj interactions .
after 
 the c36 ff was released , the namd program added the forced - based switching function , so that results 
 from c36 ff simulations in namd and charmm become consistent within 
 the differences in the simulation integrators .
the gromacs program has all the functionality needed to match with 
 the c36 ff , and its single - core code is faster than charmm and namd , 
 especially with mixed - point precision ( which used to be called single - point 
 precision ) .
the amber program is another 
 popular md code , but the latest version at the time of writing ( version 
 14 ) lacks the shift and switch cutoff methods , and thus one can only 
 use a hard truncation method for the lj term .
the force - based switching 
 function has been implemented in the development version and should 
 improve agreement with the other programs when the next version is 
 released .
openmm has been recently developed 
 as a toolkit for running md simulation with graphical processing unit 
 ( gpu ) accelerations .
however , use of the force - based switching function 
 is not straightforward , as the current nonbonded force calculation 
 in openmm only supports the potential - based switching function .
the focus of this work is to prescribe the appropriate simulation 
 protocols to run the c36 lipid ff ( and the remainder of the c36 ff 
 including the proteins , nucleic acids , carbohydrates , and 
 small organic molecules ) in a wide range 
 of available molecular simulation software packages .
since charmm , 
 namd , gromacs , amber , and openmm use different algorithms for their 
 thermostat / barostat and nonbonded interaction cutoff options , a careful 
 study on how these differences affect bilayer properties is needed . 

the importance of this effort is emphasized by a recent study in which 
 it was reported that use of the c36 lipid ff in the context of amber 
 gpus gave incorrect surface area per lipid for dppc bilayers , which 
 was due to the use of lj truncation methods that were not consistent 
 with the original parametrization . in 
 the present work ,
we systematically identify simulation protocols 
 for each program to best agree with c36 lipid ff simulations in charmm 
 and namd in terms of the following three metrics : the area per lipid 
 ( lateral density ) , the acyl chain order parameters ( lipid packing ) , 
 and the area compressibility ( dynamic lipid fluctuation ) .
the charmm - gui membrane builder(19,4042 ) now incorporates the option to generate the minimization , equilibration , 
 and production inputs ( with the optimal simulation parameters resulting 
 from this study ) for each program to allow for more general use of 
 the c36 ff . in the following sections , we describe the setup of the 
 charmm - gui input generator and the optimized protocol 
 for each program in detail . to illustrate that the optimized parameters 
 are robust and well suited for each program
, we provide the results 
 of pure lipid bilayer simulations with various lipid types ( figure 1 ) .
to determine a set of optimal simulation protocols 
 for namd , gromacs , 
 amber , openmm , and charmm / openmm simulations using the charmm c36 
 ff , a pure dppc bilayer system was built and simulated with the various 
 simulation parameters available in each program .
in addition , 1,2-dioleoyl - sn - phosphatidylcholine ( dopc ) , 1-palmitoyl-2-oleoyl - sn - phosphatidylcholine ( popc ) , 1-palmitoyl-2-oleoyl - sn - phosphatidylethanolamine ( pope ) , 1-palmitoyl-2-oleoyl - sn - phosphatidylserine ( pops ) , and palmitoylsphingomyelin 
 ( psm ) bilayer simulations were performed to validate the optimized 
 protocols derived from the dppc simulations .
the procedures of building 
 the lipid bilayer systems and the tested simulation parameters for 
 each program are provided in detail below .
all pure lipid bilayer systems were built 
 using the charmm - gui membrane builder .
a total of 80 lipid molecules were placed in each 
 lipid bilayer ( i.e. , 40 lipids in each leaflet ) with its center at z = 0 .
a water layer of 20- thickness was added above 
 and below the lipid bilayer for the dppc , dopc , popc , and pope systems 
 and that of 40- thickness for the pops and psm systems .
the 
 corresponding hydration numbers of each membrane are 34 ( dppc ) , 34 
 ( dopc ) , 32 ( popc ) , 29 ( pope ) , 67 ( pops ) , and 58 ( psm ) waters per lipid . 

kcl ions corresponding to 0.15 m were added in the pops ( including 
 neutralizing ions ) and psm systems , and no ion was added in the dppc , 
 dopc , popc , and pope systems .
while charmm and namd use the origin 
 for the periodic boundary conditions ( pbc ) , gromacs , amber , and openmm 
 use the pbc center at ( lx/2 , ly/2 , lz/2 ) where lx , ly , and lz are the system 
 lengths along the xyz directions .
therefore , the 
 whole system is translated by lx/2 , ly/2 , and lz/2 for the gromacs , amber , 
 and openmm input generation in charmm - gui .
all simulations used the c36 ff 
 for lipids and the charmm tip3p water model . to get better sampling and
check the convergence , five independent 
 md simulations were performed for each bilayer system using namd , 
 gromacs , amber , and openmm .
the simulation temperature was maintained 
 above the transition temperature of each bilayer : 300.0 ( pops ) , 303.15 
 ( dopc / popc ) , 310.0 ( pope ) , and 323.15 k ( dppc / psm ) .
pbc were employed for all simulations , 
 and the particle mesh ewald ( pme ) method was used for long - range electrostatic interactions .
the simulation 
 time step was set to 2 fs in conjunction with the shake algorithm to constrain the covalent bonds involving hydrogen 
 atoms for all programs except gromacs in which the lincs algorithm was used .
after the standard membrane 
 builder minimization and equilibration steps , the production 
 run of each simulation was performed for 250 ns .
the optimal parameters 
 were determined using the most recent version of each program ( namd 
 2.9 , gromacs 5.0 , amber14 , and openmm 6.2 ) , such that the use of previous 
 versions can cause some problems .
for example , the semi - isotropic 
 pressure coupling method was not implemented until version 6.2 of 
 openmm .
the individual simulation protocols that we tested for each 
 md program are summarized in table 1 and described in detail below .
the namd simulations were performed as a reference 
 to compare the lipid properties with those from other simulation programs . 

the force - based switching function was used for the lj interactions , 
 and three different switching ranges ( 810 , 812 , and 
 1012
hoover 
 langevin - piston method with a piston period of 50 fs 
 and a piston decay of 25 fs as well as langevin temperature coupling 
 with a friction coefficient of 1 ps were used 
 to control the pressure and temperature , respectively .
a python program was developed and used in charmm - gui 
 to convert the charmm parameter files to the namd - readable format , 
 and it is freely available in the downloadable tarball (  download.tgz  ) 
 generated by charmm - gui .
the force- and potential - based 
 switching functions 
 for the lj interactions were tested with three different switching 
 ranges ( 810 , 812 , and 1012  ) . to maintain 
 the temperature ,
a nos  hoover temperature coupling 
 method with a tau - t of 1 ps was used , and for pressure 
 coupling , a semi - isotropic parrinello 
rahman method with a tau - p of 5 ps and a compressibility of 4.5  10 bar was used . because direct 
 reading of the charmm topology and parameter files is not available 
 in gromacs , a format conversion python program was developed and used 
 in charmm - gui to generate the corresponding gromacs topology and parameter 
 files ( top and itp ) using the charmm topology and parameter files 
 ( prm , rtf , and str ) as well as the psf file .
because 
 amber does not support any switching 
 or shifting method for lj interactions , we examined the hard - truncation 
 method with three different cutoff values ( 8 , 10 , and 12  ) . 

the temperature was maintained using langevin dynamics with a friction 
 coefficient of 1 ps .
in addition , tau - p values 
 of 0.5 and 1 ps for the semi - isotropic berendsen pressure control were tested .
openmm 
 provides a library for md simulation 
 with gpu acceleration , so it does not have an executable program to 
 run md simulations .
we developed a set of python scripts that enable 
 one to run openmm simulations with simple and generalized input files , 
 and all these scripts are provided by charmm - gui .
although the standard 
 nonbonded force calculation in openmm only supports the potential - based 
 switching function , we provide an implementation of the force - based 
 switching function using openmm s customnonbondedforce class in the scripts generated by charmm - gui .
the force - based switching 
 method with the switching range of 1012  for lj interactions 
 was used for all openmm simulations ; other lj interaction parameters 
 for the dppc bilayers were tested using charmm / openmm .
langevin dynamics 
 was used for the temperature coupling , and friction coefficient values 
 of 1 and 10 ps were tested .
a semi - isotropic monte 
 carlo ( mc ) barostat method was used for the pressure 
 coupling , and pressure coupling frequency values of 5 , 10 , 15 , and 
 100 steps were examined .
openmm is able to directly read the charmm 
 topology and parameter files , so ff conversion was not necessary .
because charmm / openmm uses the openmm 
 library to calculate energies , the results between both programs should 
 be identical .
therefore , the dppc bilayer simulations to determine 
 the lj interaction parameters for both openmm and charmm / openmm were 
 performed using charmm / openmm .
the force- and potential - based switching 
 methods were tested with switching ranges of 810 , 812 , 
 and 1012  .
langevin dynamics with the friction coefficient 
 of 10 ps and the semi - isotropic mc barostat with 
 the pressure coupling frequency of 15 were used .
charmm / openmm simulations 
 for other lipids were not performed because it is identical to openmm , 
 and charmm - gui provides the charmm / openmm simulation parameters that 
 are identical to openmm simulations .
the membrane 
 properties , such as an average 
 surface area per lipid , compressibility , and deuterium order parameters , 
 were calculated and compared with available experimental values .
the 
 first 50 ns of each trajectory was discarded as an equilibration period , 
 and the properties were calculated over the last 200 ns period .
most 
 systems were equilibrated within 50 ns except for those showing the 
 phase transition ( supporting information figure s1 ) ; for such systems , the data before the phase transition 
 were discarded .
the average surface area 
 per lipid ( al ) was calculated 
 by simply dividing the system area ( axy ) by the number of lipids in each leaflet ( nl = 40 ) .
the compressibility , ka , was calculated using the 
 following equation:1where kb is the 
 boltzmann constant , t is the system temperature , 
 and al is the average of the squared fluctuation of al .
the deuterium order parameters 
 and electron density profiles ( for pops ) were calculated using st - analyzer .
the default distance cutoff of 2.4  between 
 the donor h atom and the acceptor atom was used except the lipid 
ion 
 interactions calculated with the distance cutoff of 3.8  instead . 
 an angular cutoff was not used .
the 
 bilayer properties from the dppc simulations with various simulation 
 parameters ( table 1 ) available in namd , gromacs , amber , openmm , and charmm / openmm are 
 first presented to determine the optimal simulation protocol for each 
 program using the c36 ff .
the al values from the namd simulations are 64.8 
  0.1 , 62.9  
 0.1 , and 61.6  0.1  for 810 , 812 , 
 and 1012  of the force - based switching ranges , respectively . 

these results are in good agreement with the experimental al values of 63.0  1.0 
  at 323 k. the ka values are 222.5  9.5 
 ( 810  ) , 227.1  7.1 ( 812  ) , and 242.4 
  20.9 dyn / cm ( 1012  ) , which are all in the range 
 of its experimental value of 234 dyn / cm .
although the 812  switching range appears to yield 
 the best al , the 1012 
  switching range is considered to be a standard to compare the 
 lipid properties with other programs , as 1012  is consistent 
 with the force - based switching range used to develop the remainder 
 of the additive charmm ff ( for proteins , nucleic acids , carbohydrates , 
 and small molecules ) .
we note that use of the 812  
 switching range with the full additive c36 ff is likely to be satisfactory , 
 although rigorous testing has not yet been performed using this switching 
 scheme .
we first tested the influence of mixed versus 
 double precision on the lipid properties ( supporting information table s1 ) and found no reduction in accuracy 
 with the use of the more efficient mixed precision .
second , the influence 
 of switching functions ( potential- vs force - based ) on lipid properties 
 was investigated .
al values are greatly reduced with the potential - based switching function , 
 and the 1012  switching range results in a phase transition 
 to a gel - like state . as shown in figure 2a ,
the gromacs simulation with the 1012 
  potential - based switching function ( gm - sw-3 ) shows two peaks 
 at around 49  and 58  , which represent 
 gel- and liquid - like phase states , respectively . on the other hand , 

the simulations with the force - based switching function show good 
 agreement with namd simulation results with the slightly reduced al values : 63.6  0.1 ( 810 
  ) , 61.2  0.1 ( 812  ) , and 60.3  0.1 
  ( 1012  ) .
most ka results are acceptable with the force - based 
 switching simulations showing good agreement with the namd results : 
 212.6  6.2 ( 810  ) , 213.9  7.6 ( 812 
  ) , and 212.4  6.8 dyn / cm ( 1012  ) .
the deuterium 
 order parameter values in figure 2c , d strongly suggest that the potential - based switching 
 function should not be used for the gromacs simulations with the c36 
 ff due to its disagreement with the namd results for both sn-1 and sn-2 chains . in general , these 
 gromacs simulations show higher order parameters than other simulations . 

overall , the 812  force - based switching function shows 
 best agreement with namd dppc bilayer simulations with 1012 
  force - based switching , but to be consistent with the remainder 
 of the additive charmm ff , the 1012  switching range 
 was used for the gromacs simulations of the other bilayer systems .
dppc bilayer 
 properties from some representative simulations : ( a ) 
 the distributions of surface area per lipid over the 200 ns of trajectories . 

c and 
 d ) the deuterium order parameters of ( c ) sn-1 and 
 ( d ) sn-2 chains . determining an optimal protocol for amber was 
 challenging due to the absence of any switching methods for lj interactions 
 in the current version .
the 8- cutoff shows the best al ( 62.6  0.1  ) , while the 12- cutoff shows a phase transition to 
 a gel - like state ( data is not shown ) .
this is consistent with what 
 was observed previously by gould et al . and demonstrates the need to carefully consider the cutoff scheme 
 in different programs .
compared to the namd results , the ka values are mostly in the acceptable 
 range : 235.9  12.2 ( 8  ) and 228.6  18.3 dyn / cm ( 10 
  ) .
amber simulations result in slightly reduced order parameters 
 compared to namd , but they are in the acceptable range .
there is no 
 considerable difference in between 0.5 and 1.0 ps for the pressure 
 coupling time ( tau - p ) .
thus , the 8- cutoff with tau - p of 0.5 
 ps was used for the other lipid simulations .
the cutoff methods for lj interactions 
 for both openmm and charmm / openmm were tested using charmm / openmm . 

like gromacs , the potential - based switching function reduces the al values : 61.4  0.3 ( 810 
  ) , 60.0  0.2 ( 812  ) , and 58.9  0.3 
  ( 1012  ) .
the equilibrium properties 
 of the force - based switch function systems show good agreement with 
 the namd results ; the al values are 64.6  0.2 ( 810  ) , 62.5  0.1 
 ( 812  ) , and 61.1  0.2  ( 1012 
  ) , and the ka values 
 are 227.8  14.5 ( 810  ) , 238.3  13.3 ( 812 
  ) , and 240.0  15.9 dyn / cm ( 1012  )
. the deuterium 
 order parameters of charmm / openmm simulations also show good agreement 
 with the namd results except that the potential - based switch systems 
 have higher order parameters ( figure 2c , d ) . for further openmm simulations ,
the lipid properties of the openmm simulations 
 are similar to those of the namd simulation results ( figure 2 , supporting information table s1 ) .
there is no significant difference 
 with different pressure coupling frequencies and langevin dynamics 
 friction coefficients . to reduce the gpu time ,
the pressure coupling 
 frequency of 100 steps and the langevin friction coefficient of 1 
 ps were used for other lipid simulations . the optimal protocol for each program determined from the dppc 
 bilayer simulations
is summarized in table 2 . to validate these protocols , we performed 
 the dopc , popc , pope , pops , and psm bilayer simulations and compared 
 the lipid properties of each program .
the al values are 68.7  0.1 ( namd ) , 68.0  0.1 ( gromacs ) , 
 69.7  0.1 ( amber ) , and 68.6  0.2  ( openmm ) , 
 and they are in good agreement with the experimental estimate from 
 the fully hydrated bilayer , i.e. , 67.4  1.0  at 303 k. in addition , the calculated ka values are similar among 
 the programs : 278.0  13.4 ( namd ) , 255.5  14.8 ( gromacs ) , 
 301.2  12.3 ( amber ) , and 259.5  26.9 dyn / cm ( openmm ) . 

while amber and openmm order parameters show good agreement with namd , 
 gromacs results in slightly increased order parameters for both sn-1 and sn-2 chains but they are in the 
 acceptable range ( supporting information figure s2 ) .
the experimental al of a popc bilayer is 68.3  1.5 
  at 303 k , and namd 
 shows a somewhat 
 reduced al of 65.0  
 0.1  .
gromacs results in a slightly reduced al ( 64.1  0.1  ) compared to namd , while amber shows marginally increased al ( 66.0  0.1  ) .
openmm shows the best agreement with namd in the al value ( 64.7  0.1  ) .
the ka value 
 of each program is also in the range of the namd result ( 249.6  
 19.7 dyn / cm ) : 262.5  18.2 ( gromacs ) , 270.3  18.3 ( amber ) , 
 and 274.6  31.5 dyn / cm ( openmm ) . like dopc simulation results , 

the gromacs order parameters are insignificantly increased in both sn-1 and sn-2 chains compared to other 
 programs ( supporting information figure s2 ) .
the average al value from the namd pope simulation is 57.8  0.1  , and it shows good agreement with the experimental value 
 of 59.7560.75  at 308313 k. other programs also show good agreement with 
 namd : 56.7  0.1 ( gromacs ) , 58.1  0.1 ( amber ) , and 57.7 
  0.1  ( openmm ) .
the ka values are 270.0  15.8 ( namd ) , 270.1 
  8.8 ( gromacs ) , 254.8  4.4 ( amber ) , and 266.4  24.5 
 dyn / cm ( openmm ) .
however , gromacs again yields negligibly increased 
 order parameters compared to namd , amber , and openmm ( supporting information figure s2 ) . unlike the 
 other lipid systems described above , 
 the al values of the 
 pops bilayer systems are underestimated by 10% compared to 
 its experimental al of 
 62.7  at 298 k. namd 

yields the al of 56.9 
  0.2  , and other programs show similar al values as namd : 56.0  
 0.4 ( gromacs ) and 56.9  0.5  ( openmm ) except 
 amber ( 59.0  0.1  ) .
the pops namd simulation 
 in different conditions resulted in the al of 59.9  . in this simulation , 

the system contained 100 pops lipid molecules and around 5300 water 
 molecules , and 100 k counterions were added .
thus , the slightly lower temperature with the 
 higher ion concentration in the current simulations might account 
 for the lower al values . 

the ka values are similar 
 among the programs : 261.1  22.8 ( namd ) , 262.1  33.2 ( gromacs ) , 
 270.5  28.5 ( amber ) , and 274.0  28.9 dyn / cm ( openmm ) . 
 while gromacs shows slightly increased order parameters compared to 
 namd , amber yields marginally reduced values ( supporting information figure s2 ) .
there is no significant 
 difference in the density profiles of pops bilayer systems ( figure 3 ) except that amber shows a slightly reduced hydrophobic thickness 
 ( 29.8  0.1  ) compared to namd ( 31.2  0.1  ) . 

the pops headgroup has a peak at 20  along the z - axis , and k ion shows the highest peak near 
 the pops headgroup .
( b ) the comparison of k ( solid lines ) and 
 cl ( dashed lines ) distributions among the programs . because there 
 is no reliable experimental al for the psm bilayer , the 
 gromacs , amber , and openmm al values
are compared only with the namd value . while the al values from gromacs ( 54.3 
  0.2  ) and openmm ( 55.7  0.2  ) show good agreement with the namd al value ( 55.2  0.2  ) , amber yields a 1.4  higher value ( 56.6  
 0.3  ) .
the ka values of gromacs ( 445.4  13.3 dyn / cm ) , amber ( 392.3 
  49.9 dyn / cm ) , and openmm ( 403.2  49.9 dyn / cm ) are in 
 the range of the namd ka value ( 456.4  64.8 dyn / cm ) .
the gromacs simulations result 
 in slightly increased order parameters than namd in both sphingosine and n - linked acyl chains , while amber yields slightly 
 reduced order parameters .
generally , most programs result in good agreement in the 
 order parameters , which are shown with the nmr data in figure 4 .
table 3 lists principal 
 h - bonds in the psm bilayer system and their probabilities obtained 
 from each program . in general , the h - bond probabilities are similar 
 between programs .
value sorted scd order parameters 
 of 
 psm n - linked acyl chain from nmr and each program .
the use of the c36 lipid 
 ff in various programs ( namd , gromacs , 
 amber , openmm , and charmm / openmm ) has been shown to result in comparable 
 lipid properties for dppc , dopc , popc , pope , pops , and psm bilayers 
 when the optimized simulation protocols for treatment of the nonbonded 
 interactions are applied . since each program commonly uses different 
 methods in dealing with lj cutoffs , care must be taken to choose cutoff 
 schemes that best match the approach taken during development of the 
 c36 lipid ff using the charmm program . excluding amber
, all simulation 
 programs have force - based switching functions that are required to 
 properly reproduce the published c36 ff results , and , thus , the recommended 
 cutoff scheme for these programs utilizes this approach .
the lack 
 of a switching function in amber has misled some researchers attempting 
 to apply the c36 ff in this program , although 
 the inclusion of the force - based switching function in upcoming versions 
 of amber should ameliorate these issues in the future . as we demonstrate , 
 using potential truncation at cutoffs > 8  results in gel - like 
 lipid bilayers even above the gel transition temperature .
lipid interactions will potentially result 
 in a phase change . as described in the results and 
 discussion section , the use of 8  hard truncation in 
 amber best represents the spirit of the c36 ff , although there are 
 slight discrepancies when applying this cutoff scheme to the other 
 lipids tested ( supporting information figure s2 and table s3 ) .
the discrepancies become larger when the system 
 contains ions ( pops ) or strongly interacting lipids ( psm ) . therefore
, 
 the use of the c36 lipid ff in amber ( versions 14 and earlier ) may 
 not agree with charmm / namd to the degree obtained with programs that 
 contain the force - based switching function for the lj potential .
for programs other than amber , the agreement with charmm / namd results 
 for all lipid bilayers tested in this study is satisfactory .
however , 
 there are some minor differences in that the gromacs al values are slightly reduced and thus 
 the order parameters are slightly higher .
although the cutoff scheme 
 is the same between these programs , differences in other aspects of 
 the simulation protocols , i.e. , integration scheme and temperature-/pressure - coupling , 
 can result in slight changes in lipid membrane properties .
the present 
 results indicate that the effect of different integrators and/or coupling 
 schemes on lipid bilayer properties is smaller than that due to different 
 cutoff schemes .
the usage of different cutoff schemes can also lead 
 to a bilayer phase change . while the current studies are based 
 on membrane simulations
, the 
 optimized simulation protocols should be applicable to the remainder 
 of the c36 ff including proteins , nucleic acids , carbohydrates , and 
 small molecules .
membrane simulations are especially sensitive to 
 slight changes in the lj cutoff methods , whereas other biological 
 molecules do not have the ordered self - assembly as lipids , making 
 them less sensitive to subtle details of the treatment of the lj potential . 

the optimal simulation protocol for each program has been implemented 
 in charmm - gui .
as charmm - gui will convert any patch used during the 
 system building process into appropriate simulation files , including 
 those for gromacs , the study of complex heterogeneous systems with 
 the c36 ff will be readily accessible to a range of simulation packages . 
 for example
, adding carbohydrates to proteins can be easily done with 
 charmm - gui and allows for building of such systems and their proper 
 conversion for use in gromacs .
therefore , the simulation systems and 
 inputs being seamlessly generated by charmm - gui can be used for general 
 research use in these simulation programs .
 proper treatment of nonbonded interactions 
 is essential for the 
 accuracy of molecular dynamics ( md ) simulations , especially in studies 
 of lipid bilayers . 
 the use of the charmm36 force field ( c36 ff ) in 
 different md simulation programs can result in disagreements with 
 published simulations performed with charmm due to differences in 
 the protocols used to treat the long - range and 1 - 4 nonbonded interactions . 
 in this study , we systematically test the use of the c36 lipid ff 
 in namd , gromacs , amber , openmm , and charmm / openmm . 
 a wide range of 
 lennard - jones ( lj ) cutoff schemes and integrator algorithms were tested 
 to find the optimal simulation protocol to best match bilayer properties 
 of six lipids with varying acyl chain saturation and head groups . 
 
 md simulations of a 1,2-dipalmitoyl - sn - phosphatidylcholine 
 ( dppc ) bilayer were used to obtain the optimal protocol for each program . 
 
 md simulations with all programs were found to reasonably match the 
 dppc bilayer properties ( surface area per lipid , chain order parameters , 
 and area compressibility modulus ) obtained using the standard protocol 
 used in charmm as well as from experiments . 
 the optimal simulation 
 protocol was then applied to the other five lipid simulations and 
 resulted in excellent agreement between results from most simulation 
 programs as well as with experimental data . 
 amber compared least favorably 
 with the expected membrane properties , which appears to be due to 
 its use of the hard - truncation in the lj potential versus a force - based 
 switching function used to smooth the lj potential as it approaches 
 the cutoff distance . 
 the optimal simulation protocol for each program 
 has been implemented in charmm - gui . 
 this protocol is expected to be 
 applicable to the remainder of the additive c36 ff including the proteins , 
 nucleic acids , carbohydrates , and small molecules . 
Introduction
Methods
Results and Discussion
Conclusions