membranes are integral 
 components of the cell , separating intracellular 
 compartments from the cytosol .
such membranes consist of a back - to - back 
 arrangement of lipid molecules , driven into a bilayer structure by 
 the hydrophobic effect , leaving the polar lipid head groups exposed 
 to water , and bringing the nonpolar lipid tail groups together .
the 
 composition of cell membranes is complex , with constituent species 
 including , but not limited to , saturated and unsaturated pc and pe 
 lipids , sphingomyelin and cholesterol , which serve as a matrix in 
 which membrane proteins may reside .
cell 
 membranes possess functions such as regulating transport in to and 
 out of the cell and modulating the activity of membrane embedded ion 
 channels and proteins . in order to probe the many 
 roles of membranes in the cell , membrane 
 structures
are studied experimentally using techniques such as x - ray 
 and neutron scattering , ir / raman , and nmr spectroscopy . to gain atomic - level resolution , however , membranes may also be 
 simulated computationally using molecular dynamics ( md ) .
the validity 
 of results obtained using md methods depends , to a large extent , on 
 the potential energy function , or force field , that is used .
membranes can be simulated using all - atom , united - atom , or coarse - grained 
 models , with the increasing simplicity of each representation allowing access 
 to larger models and longer time scales , at the expense of atomic 
 detail .
all - atom models may be preferred for bilayer simulations due 
 to their ability to reproduce nmr order parameters and easy combination 
 with all - atom protein , nucleic acid , carbohydrate , and small molecule 
 force fields .
one such md software package that includes all - atom 
 simulations is amber , which contains extensively validated 
 protein , nucleic , carbohydrate , and small molecule force fields .
the 
 amber simulation code has also been ported to nvidia gpu cards , making 
 simulation speeds in excess of 100 ns per day possible for systems 
 of 25,000 to 50,000 atoms , with systems of up to 4 million atoms possible 
 on the latest generation hardware .
although the amber force 
 field suite includes numerous different species of biological interest , 
 parameters for the simulation of lipids have traditionally been lacking . 

recently the lipid11 framework was introduced , which is a modular 
 lipid amber force field , allowing the simulation of a number of lipids 
 via the combination of different head and tail groups .
lipid11 used force field parameters predominantly 
 taken directly from the general amber force field ( gaff ) .
although previous studies found some success 
 in simulating lipid bilayers using gaff parameters , longer time 
 scale simulations required a surface tension term in order to keep 
 a bilayer in the correct phase at a given temperature . in this paper we draw inspiration from 
 previous work to update lipid11 headgroup 
 and tail group charges 
 and parameters to enable proper tensionless simulation of lipid bilayers , 
 thereby creating a modular amber lipid force field .
lennard - jones 
 and torsion parameters are revised , while charges are derived according 
 to the standard amber convention , as was implemented in lipid11 with 
 minor improvements in sampling for lipid14 . as such
, the resulting 
 parameters are expected to be compatible with other force fields in 
 the amber package .
six different lipid types are simulated for a total of 500 ns each , 
 with resulting area per lipid , volume per lipid , isothermal compressibility , 
 nmr order parameters , scattering form factors , and lipid lateral diffusion 
 finding favorable comparison to experiment , particularly for pc lipids . 

we also assess the conformation of the lipid tails by examining the 
 number of rotamers and rotamer sequences . to fully test the reproducibility 
 of the results ,
we believe that these parameters will be transferable thus allowing 
 the easy insertion of new lipid types into lipid14 owing to the modular 
 nature of the force field . the lipid14 force field will be released 
 with the upcoming ambertools14 .
it is our intention to incorporate 
 support for additional lipid types and other residues commonly found 
 in bilayers , as part of the upcoming release .
lipid14 aims to extend 
 the lipid11 framework to allow accurate , 
 tensionless simulation of numerous lipid types via a modular lipid 
 force field .
bond , angle , torsion , and lennard - jones ( lj ) parameters 
 in lipid11 were taken directly from the general amber force field 
 ( gaff ) .
although bond and angle parameters 
 should not require updating , it was expected that torsion and lj parameters 
 would need modification to realize this aim .
previous work on 
 lipid simulation with amber indicated that the lennard - jones and torsion 
 parameters of the lipid aliphatic tails are intricately linked to 
 the ability of the force field to reproduce experimentally observed 
 bilayer properties .
indeed , the majority 
 of all - atom lipid force field parameterization has involved the modification 
 of the lj and torsion parameters used to model the aliphatic tail 
 regions of the lipids . simulating a box of 144 pentadecane 
 chains using the standard gaff 
 lj and torsion parameters at constant pressure ( 1 atm ) and temperature 
 ( 298.15 k ) causes the hydrocarbon chains to 
freeze 
 into a crystalline state in under 2 ns ( see figure 1 ) .
such a scenario has previously been observed by klauda 
 et al . using a box of tetradecane molecules and the amber99 force 
 field .
consequently , the calculated density 
 and heat of vaporization are much higher and the diffusion much lower 
 than experiment .
a box of 144 pentadecane molecules simulated in the npt 
 ensemble 
 at 298.15 k using the general amber force field to model the carbon chains . in light of this ,
lj and torsion parameters were modified 
 to reproduce 
 the experimental density (  ) and heat of vaporization 
 ( hvap ) of alkanes covering a range 
 of chain lengths .
given that both the torsion and lj parameters affected 
 the simulated  and hvap , these parameters were altered simultaneously , with 
 the ch2ch2ch2ch2 torsion being fitted to ab initio data using paramfit .
satisfactory agreement was found by modifying only the lj parameters 
 of the hydrogen atom on the alkane chains .
the modified lj and 
 torsion parameters were tested by examining 
 the density  , heat of vaporization hvap , the diffusion d , the c nmr t1 relaxation times , and 
 the trans / gauche conformer populations 
 of a selection of hydrocarbon chains . the parameters of the 
 ester linkage region connecting headgroup 
 and tail residues in the lipids were then examined using methyl acetate 
 as a model compound .
the density and heat of vaporization of methyl 
 acetate were calculated and hvap found to be in poor agreement with experiment .
hence the lennard - jones 
 parameters of a number of atoms in this region were also adjusted 
 until better agreement with experiment was achieved .
once optimal 
 hydrocarbon and glycerol parameters were identified , 
 attention turned to the lipid partial charges .
the lipid11 force field 
 used a capping procedure , separating the lipid head and tail groups 
 into 
the standard amber resp protocol was used to generate partial charges from quantum 
 mechanical ( qm ) optimized structures , using six different orientations 
 of a single conformation .
however , in line with other 
 all - atom lipid force field parameterizations , a greater number of conformations were used per residue ( twenty - five 
 headgroup conformations , fifty tail conformations ) , with the partial 
 charges calculated as an average over all conformations .
the head 
 and tail group starting structures were extracted from previous in - house 
 bilayer simulations .
the electrostatic potential ( esp ) was calculated directly from the 
 conformations extracted from a bilayer simulation , with no qm optimization 
 performed on those structures .
charges were derived using the standard 
 amber resp protocol ( at the hf/6 - 31 g * level in gas phase ) .
finally , on identification of appropriate 
 lj parameters and calculation 
 of lipid partial charges , torsions involving the ester linkage in 
 the glycerol region were fitted to qm scans performed on a capped 
 lauroyl ( la ) tail moiety .
the 
 alkane ch2ch2ch2ch2 torsion potential was refitted using torsion scans performed 
 on 
 hexane and octane molecules using an estimation of the ccsd(t)/cc - pvqz 
 level of theory via the hm - ie relationwhere the small basis set ( sbs ) was cc - pvdz 
 and the large basis set ( lbs ) was cc - pvqz .
consequently molecules 
 were optimized at the mp2/cc - pvdz level before single - point energy 
 calculations were performed at the ccsd(t)/cc - pvdz and mp2/cc - pvqz 
 levels . in order to obtain a representative torsion fit ,
it was ensured 
 that torsion scans , conducted at increments of 15 , included 
 local minima of hexane and octane .
the paramfit program of ambertools13 was used to perform a multiple molecule weighted torsion fit , with 
 the tgt local minima of hexane and tgttt local minima of octane given a weighting of 10 , all other local 
 minima given a weighting of 4 and cis conformers 
 given a weighting of 0.1 ; all other structures were assigned a weighting 
 of 1 .
lennard - jones parameters were 
 modified while simultaneously refitting the torsion parameters ( see 
 table 1 ) until good agreement with experiment 
 was achieved for heats of vaporization and densities for a range of 
 alkane chains .
these properties were monitored by performing liquid 
 phase simulations of pentane , hexane , heptane , octane , decane , tridecane , 
 and pentadecane .
the alkane trajectories also enabled the calculation 
 of a number of other properties for comparison to experiment .
initial charges 
 for each hydrocarbon chain were generated using 
 the standard amber resp protocol ( optimization and calculation of 
 the esp at hf/6 - 31 g * level in gas phase ) .
a box of 288 ( pentane , hexane 
 and heptane ) or 144 ( octane , decane , tridecane , and pentadecane ) molecules 
 were then simulated in the liquid phase at 298.15 k for 10 ns using 
 updated lj and torsion parameters . from each liquid phase simulation
, 
 fifty different chains were extracted and used for the charge calculation , 
 with the esp calculated directly from each structure at the hf/6 - 31 g * 
 level , and partial charges derived using the resp fitting procedure . 

the 
 heat of vaporization was calculated according to1where epot(g ) is the average molecular potential 
 energy in the gas phase , epot(l ) is the average molecular potential energy in the liquid 
 phase , r is the gas constant , and t is temperature . in order to calculate the heat of vaporization , 
 two simulations were performed in a similar manner to previous work . a gas phase simulation consisting of a single 
 alkane molecule
was run for 10 ns equilibration and 50 ns production 
 under the npt ensemble , at a temperature of 298.15 k to obtain epot(g ) . a 
 liquid phase simulation consisting of a box of either 144 or 288 chains , 
 run under the npt ensemble with periodic boundary conditions using 
 particle
mesh ewald to treat long - range electrostatics and a real space cutoff of 10  , was also 
 performed for 60 ns with the first 10 ns removed for equilibration 
 to obtain epot(l ) .
the 
 temperature was maintained at 298.15 k using langevin dynamics and 
 a collision frequency of 5 ps .
pressure regulation was achieved with isotropic position 
 scaling , a berendsen barostat , and a 
 pressure relaxation time of 1 ps .
the system was heated from 0 to 
 298.15 k over 20 ps , with a force constant of 20 kcal / mol /  restraining the chains .
this restraint was gradually decreased 
 to 10 , 5 and finally 1 kcal / mol /  and the system 
 simulated for 20 ps at each value of the force constant .
the heat 
 of vaporization was then calculated using eq 1 , with results reported as block averages  standard deviation 
 using five equal blocks of 10 ns . by reducing the van der waals 
 radius ( r ) and well - depth 
 (  ) of the alkane hydrogen atom type while 
 simultaneously correcting the torsion fit , satisfactory agreement 
 with experiment
three torsion scans about the c = c double bond were also 
 performed on a cis-5-decene molecule using the mp2:cc 
 extrapolation method 
when fitting the c = c double bond torsions , the lj parameters 
 of the alkene hydrogen atom ( lipid14 atom type hb ) were scaled until 
 satisfactory agreement with experiment was found for  and hvap of cis-2-hexene , cis-5-decene , and cis-7-pentadecene .
reducing both the van der waals radius ( r ) from 1.459 to 1.25 and the well - depth (  ) from 0.015 to 0.007 resulted in far better agreement for the density 
 values (  ) ; however , it proved difficult to 
 correct the heat of vaporization ( hvap ) of cis-5-decene to that of the experimental value 
 by only modifying the lj parameters .
a similar problem has previously been encountered 
 by chiu et al . and jmbeck et 
 al . the energy profile for rotating about 
 selected torsions of a cis-5-decene molecule . energy 
 evaluated using qm and the 
 hm - ie method ( filled triangle  ) , amber with standard gaff 
 parameters ( dotted line ) , and amber with lipid14 parameters ( black 
 line ) .
the uncorrected diffusion dpbc was 
 calculated for each alkane using the slope of a mean - square displacement 
 ( msd ) plot versus time for the centers of mass , averaged over the 
 trajectories of each molecule via the einstein relation2where nf is the number of dimensions ( in this instance nf = 3 ) , and r(t ) is the distance that the molecule travels 
 in time t. it is then possible to correct for the 
 system size dependence 
 of a diffusion coefficient calculated under periodic boundaries ( dpbc ) to yield the corrected diffusion coefficient dcorr by adding the correction term derived by 
 yeh and hummer3where kb is boltzmann s 
 constant , t is the temperature ,  = 2.837297 ,  is the viscosity , and l is the length of the simulation box . the diffusion 
 was calculated from 100 ns nve ensemble simulations 
 ( extended from the 50 ns npt runs ) .
pme was used , along with a 10 
  cutoff , at a temperature of 298.15 k. in order to avoid energy 
 and temperature drift , it was necessary to remove the center of mass 
 motion every 500 steps ( nscm = 500 ) , make both the shake tolerance 
 and ewald direct sum tolerance more stringent , and reduce the time 
 step to 1 fs .
diffusion values were then calculated by taking the 
 slope of the linear 25 ns region of the msd versus time curve 
 and the correction calculated using experimental viscosity values . 

diffusion results are reported with standard deviations calculated 
 by block averaging , with each 100 ns run divided into five 20 ns blocks .
the trans , gauche , end gauche , double gauche , and kinked gtg+gtg conformer populations were 
 evaluated from the 50 ns npt runs by classifying torsion angles as 
 either gauche plus ( g+ ) 0120 , trans ( t ) 120240 , or gauche minus ( g ) 240360. c nmr t1 relaxation times 
 were calculated from the npt alkane simulations using the following 
 formula for dipolar relaxation from the reorientation correlation 
 functions of the ch vectors:4 this assumes motional narrowing and an effective c 
h 
 bond 
 length of 1.117  , with n specifying the number of protons bonded to the carbon 
 and  the ch vector .
t1 relaxation times were calculated from simulations of 
 alkanes of four different chain lengths .
these simulations were repeated 
 using the same npt alkane protocol at the experimentally relevant 
 temperature of 312 k and production runs extended to 200 ns to improve 
 sampling .
it was found that the heat of 
 vaporization of methyl acetate , a model compound representing the 
 ester linkage region connecting the lipid head and tail groups , was 
 not sufficiently close to experiment using gaff / lipid11 parameters . 
 in order to correct for this discrepancy with experiment ,
the lennard - jones 
 well - depths of the carbonyl oxygen ( oc ) , carbonyl carbon ( cc ) , and 
 ester oxygen ( os ) atoms were scaled until better agreement with experiment 
 was obtained ( see table 2 ) .
the ester oxygen 
 parameters were also applied to oxygen atoms in the phosphate region .
all values at 298.15 k. the density and heat of vaporization of methyl acetate 
 were calculated 
 using an identical procedure as for the hydrocarbon chains , with a 
 box of 288 methyl acetate molecules applied for the liquid phase calculation . 

simulations were run for 60 ns and the final 50 ns used for sampling . 

resp fitting was performed in exactly 
 the same manner to lipid11 , using the 
 capping groups shown in figure 3 .
however a 
 greater number of conformations were used to calculate the average 
 charges for each unit .
twenty - five phosphatidylcholine ( pc ) 
 and twenty - five phosphatidylethanolamine ( pe ) head groups were extracted 
 from previous bilayer simulations of dopc and pope ; while fifty lauroyl 
 ( la ) , myristoyl ( my ) , palmitoyl ( pa ) , and oleoyl ( ol ) tails were extracted 
 from bilayer simulations of dlpc , dmpc , dppc , and dopc , respectively . 

each headgroup was then capped with a methyl group and each tail capped 
 with an acetate moiety ( figure 3 ) , according 
 to the lipid11 charge derivation methodology . for each conformation
, the esp was calculated directly from each 
 structure at the hf/6 - 31 g * level using gaussian 09 .
resulting charges for the pc and pe headgroup residues 
 and the la , my , pa , and ol tail group residues are detailed in the supporting information .
two torsions involving the 
 ester linkage region ( see figure 4 ) were fitted 
 to qm data .
the scans were performed on a capped lauroyl ( la ) tail 
 residue , at 15 increments using the hm - ie method with gaussian 
 09 .
these were then fitted using paramfit for periodicity n = 1 to n = 5 
 using the genetic algorithm implemented in paramfit .
a capped lauroyl tail 
 group residue was used to fit the os - cc - cd - cd 
 and oc - cc - cd - cd torsions .
the energy profiles for rotating about selected torsions of a capped 
 lauroyl tail group residue .
energy evaluated using qm and the hm - ie 
 method ( filled triangle  ) , amber with standard gaff / lipid11 
 parameters ( dotted line ) , and amber with lipid14 parameters ( black 
 line ) .
torsion fits from the top are oc - cc - cd - cd and os - cc - cd - cd .
as can be seen from figure 5 paramfit brings 
 the os - cc - cd - cd and oc - cc - cd - cd torsions into substantially better 
 agreement with the qm data .
all values at 298.15 k. the results for the alkane properties calculated using 
 the updated 
 torsion and lj parameters are shown in table 3 .
the heat of vaporization values increase with alkane chain length , 
 following the experimental trend and converging with experiment as 
 the number of carbon atoms in the chain increases .
the simulated 
 densities are reproduced somewhat better than hvap , with an rms error of 2.60% when 
 compared to experiment .
the diffusion values lie close to experiment 
 and decrease with increasing chain length , following the experimental 
 trend ; however , the rms error between simulation and experiment remains 
 significant at 20.92% . as is also the case with the heat of vaporization 
 results ,
the main source of this discrepancy is the result for the 
 shorter alkane chains .
the better agreement with experiment of these 
 parameters at modeling longer hydrocarbon chains may arise from the 
 use of a high number of octane structures during torsion fitting . 

furthermore , this parameter set is intended for the simulation of 
 membrane lipids , which typically contain aliphatic tails that are 
 ten carbon atoms or greater in length .
thermodynamic properties for a selection of 
 alkenes calculated 
 using the updated lj and torsion parameters are shown in table 4 .
as previously stated , properties of unsaturated 
 chains are not as well reproduced as for alkanes due to the difficulty 
 in tuning lj parameters to achieve the experimental heat of vaporization , 
 resulting in an rms error of 13.80% for hvap when compared to experiment .
the fractions of trans , gauche , end gauche , double gauche , and 
 kinked gtg+gtg conformers 
 per molecule were computed for the selection of alkanes under study 
 ( see table 5 ) .
experimental data , estimated 
 by ftir , exists only for tridecane ; however , 
 the updated lipid14 parameters reproduce these results extremely well 
 with an overestimation of the end gauche and double gauche conformations only .
furthermore , the population of 
 gauche conformations per molecule falls close to the experimental 
 value of 35% for all chains investigated ( t / g ratio 
 1.86 ) , meaning that the overpopulation of trans conformations which drives gaff bilayer simulations into the gel 
 phase is avoided .
the quality of the hydrocarbon parameters 
 was further assessed 
 by calculating the c nmr t1 relaxation times for heptane , decane , tridecane , and pentadecane . 

similar to diffusion , this is a dynamic property ; however , it depends 
 on the rotation of the ch vector .
simulation values tend to converge 
 with experiment upon moving further from the end carbon of the alkane 
 chains .
although the result for pentadecane is slightly high , the 
 overall comparison between simulation and experiment is reasonable . calculated c nmr t1 relaxation 
 times for selected alkane chains and comparison to experiment .
bilayers were 
 constructed using 
 the charmm membrane builder gui at the 
 relevant experimental hydration level ( see table 6 ) and converted to lipid14 pdb format using the charmmlipid2amber.x 
 script . all systems used the tip3p water 
 model and had 0.15 m kcl salt concentration 
 added to the water layer , modeled using suitable amber parameters .
the full system was minimized 
 for 10000 steps , of which the first 5000 steps used the steepest descent 
 method and the remaining steps used the conjugate gradient method .
the system was then heated from 0 k to 
 100 k using langevin dynamics for 5 ps 
 at constant volume , with weak restraints on the lipid ( force constant 
 10 kcal mol  ) . following this
, the volume was allowed to change freely and the 
 temperature increased to a lipid dependent value ( see table 6 ) with a langevin collision frequency of  
 = 1.0 ps , and anisotropic berendsen regulation ( 1 atm ) with a time constant of 2 ps for 100 
 ps .
the same weak restraint of 10 kcal mol  was maintained on the lipid molecules .
constant pressure and constant temperature 
 ( npt ) runs were performed on the six bilayers using the amber 12 package .
the gpu implementation of the amber 12 code 
 ( bugfix 21 ) was used to run the simulations on nvidia gpu cards , achieving 
 approximately 30 ns per day for the 128-lipid bilayer systems .
bonds involving hydrogen were constrained 
 using the shake algorithm , allowing a 
 2 fs time step .
pme was 
 used to treat all electrostatic interactions with a real space cutoff 
 of 10  .
a long - range analytical dispersion correction was applied 
 to the energy and pressure .
all simulations were performed at constant 
 pressure of 1 atm and constant target temperature ( table 6 ) .
temperature was controlled by the langevin thermostat , with a collision frequency of  = 1.0 
 ps , as this method was identified as the most 
 suitable in previous work .
pressure was 
 regulated by the anisotropic berendsen method ( 1 atm ) with a pressure relaxation time of 1.0 ps .
the first 25 ns of 
 each run was removed for equilibration , resulting in a total of 500 
 ns of data per lipid system , an aggregate of 3 s of data .
results 
 are presented as block averages over the five repeats  standard 
 deviation .
the majority of analysis in this paper used ptraj or cpptraj 
 analysis routines . to check stability over 
 time of the lipid bilayer systems ,
additional gpu and 
 cpu validations were also performed , and in all cases the gpu results 
 were consistent with cpu results ( see the supporting 
 information ) .
despite the degree of uncertainty in obtaining accurate 
 experimental 
 values , the bilayer surface area each 
 lipid occupies , or area per lipid , is easily calculated from membrane 
 simulations and gives a quick indication of whether a bilayer is in 
 the correct phase at a given temperature . the area per lipid for each 
 system
the al for each lipid type is reported in table 7 , with all simulation values within 3% of experimental values , indicating 
 that all the bilayers are in the correct l-phase . 

the result for pope is closer to the older experimental al value of 56.6  than the more recent al value 5960  .
nevertheless , 
 the al should be but one of a number of 
 properties calculated to validate a lipid force field .
experimentally , the volume per lipid vl is more accurately measured and is thus a better comparison 
 for simulation results than al .
the volume 
 per lipid was calculated using the dimensions of the simulation box and the volume of a water molecule as determined 
 by simulating 1936 tip3p waters in the npt ensemble for 50 ns using 
 an identical procedure to the bilayer simulations at the relevant 
 temperature .
vl values for each 
 lipid are reported 
 in table 7 , which although systematically underestimated , 
 are within 5% of experimental values .
it is likely that the headgroup 
 lj parameters could be further tuned to remedy this discrepancy , as 
 the thorough reparamaterization of the hydrocarbon chains makes it 
 unlikely that the tail groups are causing this lack of agreement .
this intuition is confirmed when studying the lipid component volumes 
 calculated with the simtoexp software .
the headgroup volume of dopc was found to be 305.41  , which is below the experimental estimate of 319331  , while the hydrocarbon chain volume of 965.88  is closer to the experimental range of 972984  .
the volume breakdown 
 provided by simtoexp was used to calculate 
 the ratio r of terminal methyl to methylene volume . 

all lipid systems report a value of r = 1.852.17 , 
 within or very close to the experimental range of 1.82.1 .
isothermal area compressibility modulus , ka , was calculated from the fluctuation in the area per lipid . in general ,
ka values 
 fall close to experiment , with experimental values falling within 
 the standard deviation of dmpc , dppc , dopc and popc simulation results ; 
 however the pope value comes out high and there is a large standard 
 deviation in all values .
although the dopc value is above the published 
 experimental value of 300 mn m , a personal communication with e. evans revealed that this ka value has recently been revised upward to 
 318 mn m , closer 
 to the lipid14 simulation result .
the berendsen method was used 
 for pressure coupling , 
 given that it is the only barostat currently available in the amber 
 md package .
it has recently been shown that berendsen pressure control 
 is not ideal for simulations in which volume fluctuations are important , thus by implementing other barostats into amber 
 better ka results may potentially be achieved . 

this is a work in progress , the results of which will be shown elsewhere . 

furthermore , larger system sizes and longer simulation times could 
 also be investigated , as such changes have been shown to speed up 
 the convergence of ka values .
the membrane thickness was examined by 
 calculating dhh , the peak - to - peak distance , 
 from electron density profiles 
 of the membranes .
again , satisfactory agreement with experiment is 
 achieved for all lipids , though the pope value is a little high , indicating 
 that this system is slightly too ordered .
an alternative bilayer 
 thickness , the luzzati thickness db , was 
 calculated using the z - dimension of the simulation 
 box and the integral of the probability 
 distribution of the water density along the z-axis.db values are found to lie close to 
 experimental values , though the db thicknesses 
 for the saturated lipids are slightly underestimated . given that db is the distance between the points along the 
 membrane normal at which the water density is half of its bulk value , 
 this suggests that water is penetrating slightly too far into the 
 hydrophobic region of the bilayer , thereby lowering the value of db .
the 
 ordering of the lipid acyl chains may be determined by calculation 
 of the order parameter scd .
this quantity 
 can be directly compared to experimental scd values determined by h nmr or h 
d bonds with respect to the bilayer normal and was 
 calculated according to5where  is the angle 
 between the bilayer normal and the vector joining ci to its deuterium atom , and < > 
 represents 
 an ensemble average .
simulation nmr order parameters for the six lipid systems 
 and comparison 
 to experiment .
the carbon-2 atom of the sn-1 and sn-2 chains 
 display markedly different order parameters owing to the different 
 alignment of the acyl chains in this region .
experimentally , it has 
 been found that the scd order parameter 
 of the c - d bonds near the headgroup in the sn-1 chains are greater 
 than the sn-2 chains .
splitting of the order parameter 
 value of the sn-2 chain from the sn-1 chain is observed for the simulated 
 lipid systems .
the unsaturated chain of the dopc , popc and pope lipids 
 show a distinctive drop at the carbon-9 and -10 positions due to the cis double bond .
the scd values for the sn-1 chain of pope are a little high , indicating 
 that pope may be slightly too ordered . in agreement with jmbeck 

, the two chains of dopc show differing behavior , the sn-1 chain 
 having higher scd values 
 about the double bond than the sn-2 chain .
the conformation of the acyl chains may be examined by analyzing 
 the rotamers and rotamer sequences along the lipid tails and comparing 
 results to experimental data collected by fourier transform infrared 
 ( ftir ) spectroscopy .
ftir can determine 
 the number of trans ( t ) and gauche ( g ) conformers and sequences of t and g ( end gauche eg , double gauche gg , gtg , and kinks gtg ) .
the lipid bilayer simulations were analyzed 
 by denoting torsion angles  in the acyl chains as either t (  < 150 or  >
150 ) , g- ( 90  < 30 ) 
 or g+ ( 30 < 90 ) .
the rotamer sequence gtg correspond 
 to g+tg+ or g - tg- while the sequence ( or kinks ) gtg corresponds 
 to g+tg- or g - tg+ .
results are shown in table 8 and 
 are in 
 general satisfactorily close to available experimental values . the 
 discrepancies observed between simulation and experimental values 
 of gtg for
dlpc and dppc may result from 
 the experimental ambiguity in assigning gtg and gtg wagging modes .
these 
 results also confirm that the bilayers are in the correct phase , as 
 the gel - to - liquid phase transition is associated with an increase 
 in the number of gauche rotamers and kink rotamer 
 sequences .
however the eg and gg results for pope are not in accordance with the experimental values 
 obtained by senak et al .
using ftir , who 
 found a marked increase in eg , gg , and gtg+gtg values going 
 from dppe to dppc because of the tighter packing of the pe lipid in 
 the l phase .
the electron 
 density profiles 
 ( edp ) were calculated by assuming an electron charge equal to the 
 atomic number minus the atomic partial charge , located at the center 
 of each atom .
profiles have also been decomposed into contributions 
 from the following groups : water , choline ( chol ) , phosphate ( po4 ) , glycerol ( gly ) , carbonyl ( coo ) , methylene ( ch2 ) , unsaturated ch = ch and terminal methyls ( ch3 ) . 

these profiles , shown in figure 8 , are all 
 symmetrical , with water penetrating up to the carbonyl groups , leaving 
 the terminal methyl groups dehydrated in agreement with experimental 
 findings .
the electron 
 density profiles were then utilized for the calculation of scattering 
 form factors using the simtoexp software .
the 
 total and decomposed electron density profiles for each of 
 the six lipid bilayer systems with contributions from water , choline 
 ( chol ) , phosphate ( po4 ) , glycerol ( gly ) , carbonyl ( coo ) , 
 methylene ( ch2 ) , unsaturated ch = ch and terminal 
 methyls ( ch3 ) .
scattering data allow direct 
 comparison between lipid bilayer simulation and experiment , avoiding 
 any intermediate modeling of experimental raw data .
x - ray and neutron scattering form factors can be computed 
 by fourier transformation of simulation electron density profiles 
 and compared to experimental scattering data .
recent work determined 
 the areas per lipid ( ax and an ) at which dopc bilayer simulations best replicate the 
 experimental x - ray scattering and neutron scattering data by varying 
 the area per lipid through application of a surface tension , with 
 the ideal situation being ax = an = anpt ( bilayer 
 is run in the tensionless npt ensemble ) . in this work
we were concerned with validating the lipid14 parameters 
 for tensionless bilayer simulation only ; thus we report the x - ray 
 and neutron scattering form factors for anpt only .
simulation x - ray scattering form factors for the six lipid systems 
 ( black line ) and comparison to experiment ( cyan circles ) .
inset : simulation neutron scattering 
 form factors at 100% d2o ( black line ) , 70% d2o ( red line ) , and 50% d2o ( blue line ) and comparison to 
 experiment ( black , red , and blue circles , 
 respectively ) . it can be seen from figure 9 that there 
 is general agreement between both the x - ray and neutron scattering 
 form factors for all lipids for which there is experimental scattering 
 data available , indicating that the simulated bilayers have the correct 
 structure . in general the minima of the experimental f(q ) profiles 

the quantity db - h was computed from the membrane thickness values ( see 
 table 7 ) .
agreement with x - ray scattering data 
 is sensitive to the value of dhh , while 
 agreement with neutron scattering data is sensitive to the value of db .
therefore it has been proposed that bilayer 
 simulations should aim to replicate experimental db - h values to best achieve 
 agreement with both types of scattering data , where db - h = ( db - dhh)/2 .
the gromos united - atom lipid force field has been shown 
 to match experiment for simulation db - h results . as evidenced by table 7 , lipid14 db - h values are lower 
 than those found by experiment , though all simulation values do maintain 
 a large standard deviation .
in fact , analysis of db - h results for two other 
 all - atom lipid force fields , charmm36 and 
 slipids indicates that this quantity is difficult 
 to reproduce using all - atom models , with only the slipids popc result 
 falling close to experiment .
figure 10 plots 
 db - h values against area per lipid for charmm36 , slipids , and amber lipid14 , displaying 
 a downward trend in db - h with increasing area per lipid , in disagreement with 
 the experimental trend .
improving this discrepancy with experiment 
 for lipid14 may further improve the comparison with scattering data ; 
 however , present simulation scattering profiles are still seen to 
 be in satisfactory agreement with experiment .
plot of db - h versus 
 area per lipid al for the three all - atom 
 lipid force fields charmm36 ( squares ) , slipids ( diamonds ) , and amber 
 lipid14 ( circles ) . values shown for dlpc ( green ) , dmpc ( magenta ) , 
 dppc ( blue ) , dopc ( red ) , and popc ( orange ) . to assess the ability of the 
 lipid14 parameters to reproduce dynamic lipid properties ,
the lipid 
 lateral diffusion coefficient dxy was calculated using the einstein relation ( eq 2 ) with two degrees of freedom ( nf = 2 ) .
diffusion coefficients were computed for each lipid as a block 
 average over the five npt production runs .
the mean - square - displacement 
 ( msd ) curves were determined using window lengths spanning 20 ns and 
 averaged over different time origins separated by 200 ps .
the slope 
 of this curve yields the diffusion coefficient using eq 2 , with the linear 1020 ns region used to perform the 
 fit .
prior to the msd calculation , the lipid coordinates were corrected 
 to remove the artificial center of mass drift of each monolayer .
cell culture 
 membrane containing 
 78% pope at 305 k. results 
 are of the same order of magnitude as experimental values ; 
 although in general they are underestimated . unlike the bulk alkane 
 work
there is no correction term to account for collective motion 
 which can not be sampled using a periodic box of limited size .
accordingly , 
 the underestimation may be a result of size effects . as highlighted 
 by poger et al .
there is a widespread in experimental lipid lateral 
 diffusion values in the literature , with a range of experimental techniques 
 applied to the calculation of diffusion values .
even different groups applying the same experimental technique 
 do not necessarily yield comparable diffusion coefficients .
our calculated 
 diffusion coefficients are nonetheless found to be in good agreement 
 with other simulation values .
given that the production runs were performed in the npt ensemble 
 and that temperature regulation methods such as langevin dynamics , 
 which randomizes particle velocities , may affect dynamic properties 
 such as diffusion , the lipid lateral diffusion was also determined 
 in the microcanonical ( nve ) ensemble .
a single production run of each 
 lipid system was extended into the nve ensemble for 100 ns using the 
 same simulation settings as used for the alkane diffusion runs ( see parameterization strategy ) .
resulting time averaged 
 msd curves are shown in figure 11 , and calculated 
 diffusion coefficients are reported in table 9 .
although similar to the dxy values determined from the npt runs , the diffusion coefficients 
 from the nve runs are slightly higher , with the results for dppc and 
 dopc showing the largest differences .
this supports a recent study 
 on the effect of temperature control on dynamic properties by basconi 
 et al .
, who found that diffusion coefficients 
 calculated from simulations applying langevin dynamics approach the 
 coefficients derived from nve simulations , provided weak coupling 
 is used for the temperature regulation .
time averaged mean square 
 displacement of the center of mass of 
 the lipid molecules versus nve simulation time .
rattling of the lipid in the local solvation cage and slower , long distance diffusion in the 
 plane of the bilayer .
the two regimes are clearly observed in the 
 msd versus time curves ( figure 11 ) and are 
 also revealed by computing the diffusion coefficients using different 
 time ranges of the msd curve .
figure 12 plots 
 the diffusion coefficient dxy against the starting time for fitting the msd slope .
time 
 windows were either 100 ps long ( fit starts between 10 ps and 500 
 ps ) or 500 ps long ( fit starts between 500 ps and 20 ns).dxy values decrease smoothly and then converge at a value representing 
 the long time diffusion of lipids in the bilayer plane .
lateral diffusion 
 coefficients for the six lipid types calculated 
 using different time ranges of the mean square displacement curve 
 for the linear fit .
the 
 lipid14 force field represents a significant advancement for 
 the simulation of phospoholipids in the amber md package .
hydrocarbon 
 parameters have been refined , resulting in good reproduction of thermodynamic 
 and dynamic properties for a number of simple carbon chains , thus 
 we can be confident that the hydrocarbon region of the lipid membrane 
 is correctly represented .
head group parameters have also been updated , 
 with the final parameter set finding good agreement with experiment 
 for a range of properties , including the area per lipid , volume per 
 lipid , bilayer thickness , nmr order parameters , scattering data , and 
 lipid lateral diffusion , without applying a surface tension in the 
 simulations .
crucially , the experimental raw data that requires no 
 empirical input to derive , namely the nmr order parameters and scattering 
 data , are well reproduced .
results for pope however indicate that 
 pe lipids may require further attention , as the order parameter results 
 for pope indicate that it remains somewhat artificially ordered in 
 comparison to experiment .
results from five gpu repeats and cpu runs 
 are seen to be consistent ( these additional results are provided in 
 the supporting information ) , with a number 
 of tests performed on gpus using both different starting structures 
 and extending production runs to 250 ns .
future improvements may involve 
 further refinement of parameters in order to address the underestimation 
 of the volume per lipid and bilayer luzzati thickness values in addition 
 to pe lipid types .
although the present study only concerns 
 the validation of lipid14 
 for the simulation of three saturated and three unsaturated lipids , 
 the modular nature of the lipid14 force field allows for a number 
 of different lipids to be constructed from headgroup and tail group 
  building blocks and for the easy insertion of new 
 lipid species into the force field .
the lipid14 charge derivation 
 follows the usual amber convention , making this force field compatible 
 with other amber potentials , such as the general amber force field and the ff99sb protein force field . as such
, the interaction of other species , such 
 as small molecules or proteins , with lipid membranes can be studied 
 in amber using the lipid14 force field .
 the amber lipid force field has been 
 updated to create lipid14 , 
 allowing tensionless simulation of a number of lipid types with the 
 amber md package . 
 the modular nature of this force field allows numerous 
 combinations of head and tail groups to create different lipid types , 
 enabling the easy insertion of new lipid species . 
 the lennard - jones 
 and torsion parameters of both the head and tail groups have been 
 revised and updated partial charges calculated . 
 the force field has 
 been validated by simulating bilayers of six different lipid types 
 for a total of 0.5 s each without applying a surface tension ; 
 with favorable comparison to experiment for properties such as area 
 per lipid , volume per lipid , bilayer thickness , nmr order parameters , 
 scattering data , and lipid lateral diffusion . 
 as the derivation of 
 this force field is consistent with the amber development philosophy , 
 lipid14 is compatible with the amber protein , nucleic acid , carbohydrate , 
 and small molecule force fields . 
Introduction
Parameterization Strategy
Parameterization
Lipid Bilayer Simulation
Validation
Conclusions