since the first simulation 
 of an enzyme using molecular dynamics 
 ( md ) was reported by mccammon et al . in 
 1977 , md simulations have evolved to become important tools in rationalizing 
 the behavior of biomolecules .
the field has grown from that first 
 10-ps - long simulation of a mere 500 atoms to the point where small 
 enzymes can be simulated on the microsecond time scale and simulations containing millions of atoms can be considered routine . however , such simulations are numerically very intensive , and using 
 traditional cpu - centric hardware requires access to large - scale supercomputers 
 or well - designed clusters with expensive interconnects that are beyond 
 the reach of many research groups .
numerous attempts have been 
 made over the years to accelerate classical 
 md simulations by exploiting alternative hardware technologies .
some 
 notable examples include atoms by at&t bell laboratories , fastrun by columbia university and brookhaven 
 national laboratory , mdgrape by riken , and most recently anton by de shaw research llc .
all of these approaches have , however , failed 
 to make an impact on mainstream research because of their excessive 
 cost .
additionally , these technologies have been based on custom hardware 
 and do not form part of what would be considered a standard workstation 
 specification .
this has made it difficult to experiment with such 
 technologies , leading to a lack of sustained development or innovation 
 and ultimately their failure to mature into ubiquitous community - maintained 
 research tools .
graphics processing units ( gpus ) , on the other 
 hand , have been 
 an integral part of personal computers for decades , and a strong demand 
 from the consumer electronics industry has resulted in significant 
 sustained industrial investment in the stable , long - term development 
 of gpu technology .
in addition to low prices for gpus , this has led 
 to a continuous increase in the computational power and memory bandwidth 
 of gpus , significantly outstripping the improvements in cpus . as a 
 consequence
, high - end gpus can be considered standard equipment in 
 scientific workstations , which means that they either already exist 
 in many research laboratories or can be purchased easily with new 
 equipment .
this makes them readily available to researchers and thus 
 attractive targets for acceleration of many scientific applications 
 including md simulations . the nature of gpu hardware , however , 
 has until recently made their 
 use in general purpose computing challenging to all but those with 
 extensive three - dimensional ( 3d ) graphics programming experience . 

however , the development of application programming interfaces ( apis ) 
 targeted at general purpose scientific computing has reduced this 
 complexity substantially such thatgpus are now accepted as serious 
 tools for the economically efficient acceleration of an extensive 
 range of scientific problems .
the computational 
 complexity and fine grained parallelism of md 
 simulations of macromolecules makes them an ideal candidate for implementation 
 on gpus . indeed , as we illustrate here for implicit solvent and in 
 a subsequent paper for explicit solvent
, 
 the careful implementation of modern md algorithms on gpus can provide 
 capability , in terms of performance , that exceeds that achieveable 
 with any current cpu - based supercomputer .
several previous studies 
 have investigated the use of gpus to accelerate md simulations . for a detailed review of the use of gpus for acceleration of condensed 
 phase biomolecular md simulations
we 
 present our high - performance gpu implementation 
 of implicit solvent generalized born ( gb ) md for the amber and charmm pairwise 
 additive force fields on cuda - enabled nvidia gpus .
we have implemented 
 this within the amber pmemd dynamics engine in a manner 
 that is designed to be as transparent to the user as possible , and 
 we give an overview of what the code currently supports , as well as 
 our plans for future developments .
we discuss the specifics by which 
 we exploit the processing power of gpus , both in serial and using 
 multiple gpus , and show the performance that can be achieved in comparison 
 to conventional cpu clusters .
we also discuss our implementation and 
 validation of three specific precision models that we developed and 
 their impact on the numerical results of implicit solvent md simulations .
as illustrated 
 by figure 1 , gpus offer a 
 tremendous amount of computing power in a compact package .
this , however , 
 comes at the cost of reduced flexibility and increased programming 
 complexity as compared to cpus .
in order to develop software that 
 runs efficiently on gpus , it is necessary to have a thorough understanding 
 of the characteristics of the gpu hardware architecture .
a number 
 of manuscripts have already discussed this in detail in the context 
 of md .
for this reason , we 
 provide simply a brief overview of the complexities involved in programming 
 gpus as they relate to our implementation , focusing on nvidia hardware . 
 for a more detailed description ,
peak floating - point operations per second ( flop / s ; left ) 
 and memory 
 bandwidth ( right ) for intel cpus and 
 nvidia gpus .
a gpu is an example 
 of a massively parallel stream - processing architecture which uses 
 the single - instruction multiple data ( simd ) vector processing model . 
 unlike a regular cpu , which typically operates on one to four threads 
 in parallel , gpus typically process threads in blocks ( termed warps 
 within the cuda programming language ) 
 containing between 16 and 64 threads .
these thread blocks logically 
 map to the underlying hardware , which consists of streaming multiprocessors . 
 at the time of writing , high - end gpus typically have between 16 and 
 32 multiprocessors . for example , an nvidia m2090 gpu consists of 16 
 multiprocessors , each containing 32 cores for a total of 512 cores . 

all threads in a single block must execute the same instruction on 
 the same clock cycle .
this necessarily implies that , for optimum performance , 
 codes must be vectorized to match the size of a thread block . branching 
 must therefore be used with extreme care since if any two threads 
 in the same warp have to follow different code paths of the branch , 
 then threads in the warp will stall while each side of the branch 
 is executed sequentially .
the memory hierarchy 
 of gpus has its origins in their graphics lineage , and the high density 
 of arithmetic units comes at the expense of cache memory and control 
 units .
all of the cores making up a multiprocessor have a small number 
 of registers that they can access , a few kilobytes ( 64 kb on an m2090 ) 
 of shared memory [ this can be split into directly accessible memory 
 and l1 cache ; in the case of an m2090 , it can be split 48/16 kb or 
 16/48 kb ; in the case of amber , the configuration is switched at runtime 
 for optimal performance of a given kernel ] which is private to each 
 multiprocessor and a small amount ( typically 48 kb ) of high - speed 
 but read - only texture memory .
the majority of the memory ( 6 gb on 
 an m2090 ) , termed global device memory , is available to all multiprocessors . 

while being fast compared to the main memory accessible by cpus , access 
 to the device memory by gpus is still relatively slow compared to 
 the local cache memory .
the nature by which the multiprocessors are 
 connected to this memory also means that there is a significant performance 
 penalty for nonstride-1 access .
finally , it should be noted that currently the cpu and gpu memories are in different address spaces and this requires careful consideration . 

the unique nature of this memory model leads to several considerations 
 for optimizing gpu performance , including optimizing device memory 
 access for contiguous data , utilizing the multiprocessor shared memory 
 to store intermediate results or to reorganize data that would otherwise 
 require nonstride-1 memory accesses , and using the texture memory 
 to store read - only information , such as various force field parameters , 
 in a fashion that allows very rapid access .
as mentioned 
 above , the cpu and gpu memories are , at the time of writing , in different 
 address spaces .
this means it is up to the programmer to ensure that 
 the memories are synchronized as necessary to avoid race conditions . 

however , there is a big performance penalty for such synchronizations 
 which have to occur via the peripheral component interconnect express 
 ( pcie ) bus , and thus they should be avoided unless absolutely necessary .
the traditional 
 method for programming scientific algorithms in paralleluses the 
 message passing interface ( mpi ) in which 
 each thread runs in a separate address space .
when running gpus in 
 parallel under an mpi paradigm , additional complexity is introduced 
 since sending data between two gpus involves copying the data from 
 the memory of the sending gpu to the cpu memory of the corresponding 
 mpi thread over the pcie bus , an mpisend by this cpu thread , and corresponding 
 mpireceive by the receiving cpu thread , which copies the data between 
 the memories of the cpus , and finally copying the data to the memory 
 of the receiving gpu .
clearly , this introduces additional considerations 
 for maximizing parallel performance as compared to traditional cpu 
 programming . at the time of writing
, there are efforts to streamline 
 gpu to gpu communication , particularly within a single node but also 
 for infiniband connections between nodes .
one such approach under 
 development by nvidia and mellanox is termed gpudirect , which ultimately seeks to unify address spaces 
 between multiple cpus and gpus . currently , the degree to which this 
 can be utilized is heavily dependent on the underlying hardware design . 
 therefore , at present , the added complexity of using the advanced 
 features of gpudirect , beyond the pinned memory mpi optimizations 
 offered by gpudirect version 1 , on the large number of possible different 
 hardware combinations is not worth the effort for a widely used production 
 code .
early versions 
 of gpus in nvidia s lineup ( prior to the gt200 model ) only 
 supported single precision ( sp ) floating point arithmetic .
this was 
 due to the fact that graphics rendering did not require double precision 
 ( dp ) .
scientific algorithms , however , typically require dp arithmetic 
 ( for a discussion in the context of quantum chemistry , see for example 
 the work by kniziaet al . ) .
the generation 
 of gpus at the time of our initial implementation ( 2008 ) supported 
 dp in hardware , but only at 1/8 the performance of sp . in the latest 
 generation of cards , at the time of writing , termed the fermi lineup 
 by nvidia , the dp to sp performance ratio is 1/2 and thus equivalent 
 to that in cpus .
this , however , only holds for the professional ( termed 
 tesla ) series of cards .
the significantly cheaper gaming cards ( termed 
 geforce ) still only support dp at a fraction of the speed of sp .
it 
 is therefore important to optimize the use of dp such that it is only 
 used when necessary to maintain numerical accuracy .
early use of gpus 
 for scientific computing was hampered by the lack of an application 
 programming interface ( api ) for general purpose calculations .
the 
 problems to be solved had to be described in terms of a graphics pipeline 
 employing either opengl or directx , which made the software development 
 time - consuming and hardware - specific .
the barrier to utilizing gpu 
 hardware for general purpose computation has since been reduced by 
 the introduction of gpu programming models such as the brook stream 
 programming language , opencl , and nvidia s compute unified device architecture 
 ( cuda ) and the availability of corresponding 
 software development toolkits ( sdks ) .
the amber implementation uses 
 cuda , which is a relatively simple extension of the standard c programming 
 language that allows one to code in an inherently parallel fashion 
 and perform all necessary operations to access and manipulate data 
 on a gpu device .
realizing the full potential of gpus , however , still 
 requires considerable effort as indicated above and outlined below 
 to take advantage of the particular gpu architecture , and not all 
 algorithms are suitable to achieve good performance on these massively 
 parallel processors .
the nature of md simulations requires what in computer science 
 is referred to as strong scaling , that is , reduction of the solution 
 time with an increasing number of processors for a fixed total problem 
 size .
this enables access to simulations at longer time scales , which 
 is required for a proper convergence of results .
this becomes more 
 important as one moves to larger system sizes since the number of 
 degrees of freedom increases .
weak scaling , that is , the solution 
 time with the number of processors for a fixed problem size per processor , 
 is only of secondary importance , since this merely enables simulating 
 larger molecules at currently attainable time scales .
our implementation 
 therefore has focused on accelerating problem sizes that correspond 
 to those typically studied by amber users . in the case of gb simulations , 

the initial 
 driving force in accelerating amber implicit solvent 
 gb calculations with gpus was to provide the scientific community 
 with a computational tool that would allow an individual researcher 
 to obtain performance on a simple desktop workstation equivalent to 
 that of a small cpu cluster .
such a tool alleviates the costs , both 
 capital and recurring , involved in purchasing , maintaining , and using 
 individual research compute clusters . to this end , our goal was that 
 a single state - of - the - art gpu should provide a performance equivalent 
 to that of four to six high - end cpu cluster nodes .
such an approach 
 also removes the need to purchase and maintain expensive interconnects 
 that are required to achieve scaling even on a modest number of nodes . beyond this initial serial development , which was first released 
 as an update to amber 10 in august 2009 , 
 we have also developed a parallel implementation based on the mpi-2 message passing protocol , released as an update 
 to amber 11 in october 2010 , that allows 
 a single job to span multiple gpus .
as shown below , it is possible with this 
 implementation to achieve a performance improvement that goes beyond 
 simply making a desktop workstation faster , ultimately providing a 
 performance capability that surpasses what is achievable on all current 
 conventional supercomputers . achieving this level of performance required 
 implementing the entire implicit solvent md algorithm including energy 
 and force evaluations , restraints , constraints , thermostats , and time 
 step integration on the gpu .
as described in section 3.2 , cpu to gpu communication only occurs during i / o or to 
 some extent when data is sent between gpus during parallel runs . while we have designed our gpu implementation to achieve substantial 
 acceleration of implicit solvent md simulations over that achievable 
 with amber s cpu implementation , our overriding goal has always 
 been to maintain the precision of the calculations . to this end
, we 
 have focused on ensuring that gpu simulations will match cpu simulations . 

all approximations made in order to achieve performance on gpu hardware 
 have been rigorously tested as highlighted in the following sections .
an additional design goal has been to attempt to preserve forward 
 compatibility of our implementation .
using the cuda programming language 
 provides this by abstracting the program from the underlying hardware . 

the gpu accelerated version of amber can be used on all nvidia cards 
 that support double precision in hardware , that is , those with hardware 
 revision 1.3 or 2.0 or higher .
our choice of cuda and nvidia graphics 
 cards was largely guided by the fact that , at the time we began this 
 work , opencl was not mature enough to offer the same performance and 
 stability benefits that cuda did .
a port to opencl is certainly possible , 
 and this would support amd hardware .
however , with the public release 
 of the cuda api by nvidia and the release 
 of cuda compilers for x86 platforms by pgi , it is possible that an amber implementation will soon be available 
 on a variety of accelerator hardware .
we 
 have attempted to make our gpu implementation all - inclusive of the 
 features available in the pmemd program . at the time of writing ,
the 
 majority of features applicable to implicit solvent simulations are 
 available as described below .
support is provided for all gb models 
 currently implemented within amber as well as the analytical linearized poisson 
in addition to constant energy simulations , 
 thermostats have been implemented to perform constant temperature 
 simulations .
this includes all three thermostats available in pmemd , 
 that is , the berendsen weak coupling algorithm , the andersen temperature coupling scheme , and the langevin dynamics thermostat .
constraints for hydrogen bond distances use a gpu version 
 of the standard shake algorithm employed in pmemd , 
 and harmonic restraints to a reference structure are supported . to the best of our knowledge
, no gb formalism currently exists that 
 corrects for the errors introduced by the use of cutoffs for long - range 
 nonbonded interactions .
the use of cutoffs in gb simulations as implemented 
 in pmemd does not conserve energy , and their use involves an approximation 
 with an unknown effect on accuracy .
for this reason , we chose not 
 to implement van der waals ( vdw ) and electrostatic cutoffs in the 
 gpu version of this code .
[ cutoffs for the nonbonded interactions 
 are implemented for explicit solvent simulations with periodic boundary 
 conditions using the particle mesh ewald ( pme ) method , as described 
 in a later paper .
a design feature of the gpu code that 
 goes beyond the cpu implementation is the deterministic nature of 
 the implementation on a given hardware configuration .
serial cpu calculations 
 for a given set of input parameters on identical hardware are perfectly 
 reproducible .
this does not hold for the parallel cpu implementation 
 since the need to load balance aggressively to achieve good parallel 
 scaling means that the order of numerical operations is not defined , 
 and therefore two simulations started from identical conditions will 
 always diverge due to rounding differences .
this poses a problem when 
 transitioning to microsecond or greater simulation time scales since 
 it can be of advantage to store trajectory information less frequently 
 than what is optimal in order to conserve available storage space 
 and produce data files of manageable size .
it is thus not possible 
 to go back to a given point of the simulation and analyze the trajectory 
 in finer detail by restarting and sampling more frequently unless 
 the implementation is deterministic .
the deterministic nature of the 
 gpu code coupled with machine precision binary restart files ( currently 
 under development ) makes this mode of simulation possible .
another key feature and a primary design 
 goal of our gpu acccelerated implementation is that its use is completely 
 transparent to the user .
as far as the user is concerned , our gpu 
 implementation is indistinguishable from the cpu implementation , and 
 using the gpu version of the code is simply a case of switching the 
 executable name from pmemd to pmemd.cuda or from pmemd.mpi to pmemd.cuda.mpi for the mpi parallelized implementation .
all other items such as 
 input and output files and regression tests within the code remain 
 identical . the only difference to be noticed by the user is an increase 
 of performance .
this guarantees effective uptake of our gpu implementation 
 by the scientific commmunity because there is no learning curve for 
 the use of the code , and all tools and scripts that have been developed 
 for the cpu version of pmemd can be utilized without modifications .
the maximum system size that can be treated 
 with the gpu implementation is a function of both the gpu hardware 
 and the md simulation parameters .
in particular , langevin temperature 
 regulation and the use of larger cutoffs for the effective born radii 
 calculations increase the memory requirements .
the physical gpu hardware 
 also affects memory usage since the optimizations used are nonidentical 
 for different gpu types .
table 1 gives an overview 
 of the approximate maximum atom counts that can be treated with the 
 present version of the code .
the dominant sources of gpu memory usage 
 are the output buffers used for the nonbonded interactions as described 
 in section 3.2 .
the memory used by those 
 buffers is proportional to the square of the number of atoms
. currently , 
 the atom count limitations imposed by gpu memory usage are roughly 
 identical in serial and parallel .
all simulations use shake ( amber input ntf=2 , 
 ntc=2 ) ; a time step of 2 fs ; the hawkins , cramer , truhlar gb model ( amber input igb=1 ) ; the default cutoff value 
 of 25  for gb radii ( amber input rgbmax=25 ) ; and temperature 
 control with the langevin thermostat ( amber input ntt=3 ) , if applicable . 

error - correction code ( ecc ) was switched off on the tesla cards . in classical md
, the majority of the computational effort is spent 
 evaluating the potential energy and gradients , which has to be repeated 
 each time step . in the case of the amber pairwise
additive force fields , the potential takes the form1where the bond and angle terms are represented 
 by a simple harmonic expression with force constants bi and ai and equilibrium bond distances / angles ri , eq and i , eq , respectively .
torsional potentials for the dihedral angles 
 are represented using a truncated fourier expansion in which the individual 
 terms have a potential vi , n with periodicity n and phase shift 
 i , n .
the last 
 two terms are the vdw interaction represented by a lennard - jones potential 
 with diatomic parameters aij and bij and the 
 electrostatic interaction between atom - centered point charges qi and qj separated by the distance rij .
the prime on the summation of the 
 nonbonded interactions indicates that vdw and electrostatic interactions 
 are only calculated for atoms in different molecules or for atoms 
 in the same molecule separated by at least three bonds .
those nonbonded 
 interactions separated by exactly three bonds ( 14 interactions ) 
 are reduced by the application of independent vdw and electrostatic 
 scale factors , termed scnb and scee in amber , which are dependent 
 on the specific version of the force field ( 2.0 and 1.2 , respectively , 
 for the ff99sb version of the amber force 
 field ) .
the charmm force field takes 
 a similar form but includes three additional bonded terms:2the two - body urey 
bradley terms account 
 for angle bending between atoms that are separated by two bonds using 
 a harmonic potential with force constant ui and equilibrium distance ri , eq .
the improper dihedrals with force 
 constant ki describe 
 out - of - plane bending and are used to maintain planarity and prevent 
 undesired chiral inversions . within the amber force field ,
the third 
 additional term is a cross term between two sequential protein backbone 
 dihedral angles , termed cmap .
additionally , the charmm and amber force fields handle 14 
 nonbonded interactions in a different manner .
the single prime on 
 the electrostatic summation has the same meaning as described above 
 for the amber force field with the exception that 14 interactions 
 are not scaled .
the double prime on the vdw summation implies the 
 same exclusions as the single prime but the use of different values rijmin and ij for 14 interactions . in the gb implicit solvent model
, 
 the effect of a surrounding solvent 
 is described via a continuum electrostatics model that uses a pairwise 
 descreening approximation and in general also includes a debye  hckel 
 term to account for salt effects at low salt concentrations .
the general 
 form of the correction to the energy of the solute is given as3where qi and qj are the 
 atomic partial charges , r is the relative permittivity 
 of the solvent , and  is the debye 
the function fijgb interpolates between an effective born radius ri when the distance rij between atoms is short and rij at large distances according 
 to4the effective born radius ri reflects how deeply buried a charge 
 is in the low - dielectric medium ( solute ) , and it depends on the intrinsic 
 radius i of an atom and the relative 
 positions and intrinsic radii of all other atoms in the solute .
the 
 models implemented in pmemd differ in the intrinsic radii and the 
 function that is used to determine the effective born radii .
the remainder of this section discusses key aspects of the gpu 
 implementation of these equations .
this includes design features to 
 maximize performance while preserving accuracy as well as addressing 
 issues of reproducibility which are critical when moving to longer 
 time scale simulations .
one of the 
 initial ideas we considered when implementing gpu acceleration for 
 amber was to write an entirely new code from scratch in a combination 
 of c++ and cuda .
however , this was quickly dismissed for a number 
 of reasons . in particular , we wanted to1.keep the maintenance of the amber code 
 base simple2.minimize 
 the amount of coding required3.simplify the way in which features 
 are ported to the gpu4.maintain backward compatibility with 
 existing input files and regression tests . keep the maintenance of the amber code 
 base simple minimize 
 the amount of coding required simplify the way in which features 
 are ported to the gpu maintain backward compatibility with 
 existing input files and regression tests .
the approach we ultimately chose was to utilize the 
 existing fortran code base in pmemd and extend it with calls to specific 
 cuda kernels for the gpu acceleration with the cuda code protected 
 with # ifdef cuda preprocessor directives .
a series 
 of gpu synchronization routines provide an abstract way 
 to copy relevant data to and from the gpu memory , for example , gpu_upload_vel ( ) 
 and gpu_download_vel ( ) for atomic velocities .
a complete list of synchronization 
 routines is provided in the supporting information . for performance
, we have implemented the entire md algorithm on 
 the gpu , which means uploads ( to gpu memory ) are only needed at the 
 beginning of a run and downloads ( to cpu memory ) are only needed when 
 i / o is required , for example , downloading the coordinates to write 
 to the trajectory file .
it remains nevertheless simple to add new 
 features to the code in a way that will work for both the cpu and 
 gpu ( see the supporting information for 
 an example ) . while this limits the performance for the new feature 
 due to repeated up- and downloads , it does provide a mechanism by 
 which new features can be easily added and tested with the gpu code 
 before writing an additional gpu kernel . in order to achieve performance in 
 parallel on cpus ,
it has always been necessary to include complex 
 dynamic load balancing algorithms and extensive use of asynchronous 
 communication within the implementation .
this arises from the fact 
 that the cpus are expected to perform a number of tasks above and 
 beyond running the underlying simulation .
for example , the operating 
 systems run multiple daemons controlling i / o and various other os 
 related tasks .
in addition , the latency between any two cpus in a 
 large parallel run can span an order of magnitude or more . as a result
a downside of this is that the order of operations is not well - defined 
 in the cpu implementation .
this results in rounding differences between 
 otherwise identical runs . due to the chaotic nature of numerically 
 integrating newton s equations of motion
, this means that two 
 initially identical simulations will always decorrelate on time scales 
 of a few nanoseconds or less .
there is nothing inherently wrong with 
 this ; the two simulations are just exploring two different regions 
 of phase space .
as mentioned previously , this does , however , pose 
 a problem when one starts to routinely run simulations on the microsecond 
 time scale since it can be necessary to return to a given point in 
 a simulation and repeat it exactly while saving the output more frequently . 

. attempted to address with a bit - wise 
 reversible integrator based on scaled 
 integers .
the architectural design of the gpu with many floating 
 point units controlled by a single instruction unit , and the fact 
 that the gpu is not required to time slice with os related tasks , 
 means that work can be divided between gpu threads and indeed entire 
 gpus in a predetermined fashion without detrimentally impacting performance . 

this can be exploited by careful programming and bookkeeping to ensure 
 that the order of floating point operations is always predefined at 
 each given time step .
the gpu implementation of pmemd has been designed 
 such that any two simulations with identical starting conditions run 
 on the same gpu model will always be bit - wise identical .
this is extremely 
 useful for debugging purposes , for example , for detecting shared memory 
 race conditions , and has also been used to determine that the use 
 of ecc on gpus has little to no impact on the reliability of amber 
 simulations . for anderson and langevin 
 thermostatted simulations ,
it is necessary for the random number generator 
 ( rng ) to also be perfectly deterministic in order for any two initially 
 identical simulations to be reproducible . to this end
, we use the 
 parallelized rng that is implemented in the curand library and available 
 with the cuda toolkits . in amber ,
all of the traditional cpu 
 codes are written entirely using double precision ( dp ) floating point 
 arithmetic .
uses single precision ( sp ) floating point numbers 
 throughout the calculations with the exception of a single double - precision 
 accumulator in the reduction phase of the force accumulation .
use 
 of sp in all places within the code , however , can cause substantial 
 instabilities in the md simulations .
while this error can be 
 hidden using tightly coupled thermostats , the true effects of such 
 approximations have not been well characterized .
we distinguish 
 three different precision models in our gpu implementation in which 
 the contributions to the nonbonded forces are calculated in single 
 precision arithmetic , but bonded terms and force accumulation are 
 in double precision ( spdp ) , or everything is computed and accumulated 
 in single precision ( spsp ) or double precision ( dpdp ) . the exception 
 to this is the shake algorithm ( see below ) , which is implemented in 
 dp for all precision models since it involves calculating relative 
 differences in distance on the order of 10  . 
 attempting to use sp for the shake algorithm as implemented leads 
 to numerically unstable simulations .
the aim in developing our spdp 
 precision model for the gpu implementation of pmemd was to achieve 
 numerical stability during md simulations equivalent to that of traditional 
 dp implementations but with performance as close to the spsp model 
 as possible .
as highlighted in the subsequent sections , the spdp model 
 achieves these aims and for this reason is the default precision model 
 used in the gpu implementation of pmemd , although the other precision 
 models can be chosen if desired .
our approach to the calculation 
 of nonbonded interactions is similar to that described in friedrichs 
 et al .
having been developed at the same 
 time and with overlapping authors but with several differences and 
 additional optimizations , including the way in which accumulations 
 are handled .
the algorithm used for the calculation of the nonbond 
 forces is identical for all three precision models .
the pairwise interactions 
 between atoms i and j , which can 
 schematically be represented by a matrix as in figure 2 , are grouped together into tiles of dimension w  w. the evaluation of the nonbond forces 
 for each tile is dynamically assigned to independent warps across 
 all of the streaming multiprocessors ( sms ) of a given gpu , which are 
 each running a sufficient number of threads to bury operational latency . 

the reason for doing this is that the gpu schedules work in terms 
 of warps which effectively perform the same mathematical operation 
 on w values at once . in the case of nvidia ,
all current 
 gpus process warps as 32-threads - wide ; in this case , optimum performance 
 is achieved by making w = 32 .
since the interaction 
 between two atoms i and j is symmetric , 
 only the blue diagonal tiles and the green off - diagonal tiles need 
 to be considered for the calculation as described below .
unlike the 
 cpu code , which can just simply loop over atoms i and j , it is crucial to divide up the calculation 
 on the gpu into as many warp size blocks as possible . however , this 
 presents a problem when it comes to accumulating the forces for each 
 atom since the warps can be scheduled for computation in any order . 

a nave accumulation of the resulting forces on each 
 atom can thus lead to race conditions and incorrect results .
the use 
 of atomic operations within the nonbond routine to avoid this problem , 
 however , represents a serial bottleneck that would be unacceptable 
 in terms of performance .
the approach we use is thus to allocate ( n / w ) + 1 output buffers per atom where n is the total number of atoms and w is 
 the warp size . in the case of the nonbond force calculation ,
each 
 output buffer is three double - precision - values - wide for the spdp and 
 dpdp precision models and three single - precision - values - wide for the 
 spsp precision model corresponding to the force in the x , y , and z directions of cartesian 
 coordinate space .
forces calculated within a tile can then be summed 
 for each atom without fear of overwriting memory effectively providing 
 a separation in space between tiles . at the conclusion of the force 
 calculation
, we assign one thread per atom in linear order and cycle 
 through the output buffers , summing them to obtain the total force 
 on each atom .
schematic representation of the work - load distribution 
 for the 
 calculation of nonbond forces with n atoms .
each 
 square represents the interactions between two atoms i and j for which the resulting forces need to be 
 evaluated .
these are grouped together in tiles of size w  w that are each assigned to an independent 
 warp . due to symmetry ,
only the blue diagonal tiles and the green 
 off - diagonal tiles need to be considered for the calculation . for 
 details , see the text .
while this approach provides a separation in space 
 between tiles , 
 it is still necessary to deal with race conditions within a tile . 

we have to distinguish on - diagonal tiles ( blue ) and off - diagonal tiles 
 ( green ) , see figure 2 .
the on - diagonal tiles 
 include the interaction of atoms i to i + w  1 with themselves and thus include 
 the self - interaction from the gb term .
in contrast , off - diagonal tiles 
 include the interaction of atoms i to i + w  1 with j to j + w  1 , where i  j. to avoid race conditions , it is necessary 
 to handle these two types of tile in different ways . on -
diagonal 
 tiles store two copies of the coordinates and associated 
 parameters for each atom , one set in shared memory ( index atom i in figure 2 ) and the second set 
 in the registers ( index atom j in figure 2 ) .
each thread in the warp then runs linearly through 
 shared memory ( atom j ) , accumulating only the force 
 acting on the register atoms i in the corresponding 
 force output buffer . as illustrated in figure 2 , on the first iteration , the force on atom i resulting 
 from interacting with itself is obtained and stored in the corresponding 
 register for i , while simultaneously the force on 
 atom i + 1 resulting from interacting with atom i is obtained and stored in the register for atom i + 1 . similarly the forces on atoms i + 
 2 to i + w  1 resulting 
 from interacting with atom i are obtained and stored 
 in the corresponding output buffer .
all threads then step to the next 
 column to obtain the forces on atoms i to i + w  1 due to interactions with 
 atom i + 1 etc .
this approach does some excessive 
 computational work since , for any given atom pair ij , the force on atom j is just the negative of the 
 force on atom i , and thus only the upper or lower 
 triangle of the tile would have to be considered .
however , this is 
 far outweighed by the advantage of avoiding race conditions that would 
 result from several threads updating their force contribution to the 
 same atom . in the case of off - diagonal tiles , a different approach 
 is taken 
 such that the symmetry in the interactions is exploited while avoiding 
 race conditions .
again , the coordinates and associated parameters 
 of atoms labeled j to j + w  1 are placed in the registers , while the coordinates 
 and associated parameters of atoms i to i + w  1 are placed in shared memory .
if 
 the off - diagonal tiles were handled in the same fashion as the on - diagonal 
 tiles , on the first iteration , the force on atom j due to interacting with atom i would be calculated 
 and from this the corresponding force on atom i obtained 
 by negation . at the same time , the force on atom j + 1 due to interacting with atom i would be calculated 
 and from this the corresponding force on atom i obtained 
 by negation .
thus , the forces on atom i due to all 
 atoms j of the tile would be calculated at the same 
 time , which would require atomic operations to accumulate correctly . 

the solution to this problem that is implemented in amber is to partition 
 the off - diagonal tiles in time as illustrated in figure 2 . by starting each thread offset by its thread
i d , we avoid 
 race conditions in the accumulation and eliminate the need for performance 
 destroying atomic operations . the gb terms are calculated 
 in an identical fashion to the nonbond terms described above . the 
 born radii are first calculated within their own kernel and reduced 
 to per atom radii in an analogous fashion to how the nonbond forces 
 are handled .
the bonded and 14 
 terms represent a very small fraction of the computational workload 
 ( typically < 1% of an iteration ) , and thus optimization of their 
 calculation is not critical for performance compared with , for example , 
 the nonbond interactions
. however , efficient calculation on a gpu 
 still requires some consideration in how these terms are calculated 
 in order to exploit the massive parallelism within the gpu while avoiding 
 race conditions and memory overwrites .
our implementation simply creates 
 a list of the interactions , sorted by type , and then divides up the 
 bonded and 14 terms across sms on a per interaction basis . 

since the resulting forces need to be summed for each atom , there 
 is the potential for a race condition to occur during this reduction 
 if two or more gpu threads attempt to accumulate forces for the same 
 atom at the same time in global memory .
our initial attempts to avoid 
 this used atomic operations for reduction of the resulting forces 
 to individual atoms but showed very poor performance .
our solution 
 to this is to make use of the tile reduction buffers from the nonbond 
 calculation as described above and then sum these forces in the reduction 
 step used in the nonbond calculation . at the time of writing ,
these are handled in an identical fashion to the bonded interactions 
 being calculated as a bond between an atom and a fixed virtual particle 
 representing the reference structure .
the shake implementation is currently 
 restricted to hydrogen atoms since this represents > 99% of the 
 types 
 of simulation amber users run .
this restriction has the benefit that 
 each heavy atom and its attached hydrogen atoms can thus be handled 
 independently .
our approach is to simply spawn a thread for each heavy 
 atom during the shake calculation .
the use of double precision entirely within 
 the shake routine provides sufficient precision to handle shake tolerances 
 on the order of 10  and thus avoids complexities 
 involved in attempting to implement shake constraints in an internal 
 coordinate representation as would be required if shake was to be 
 carried out in single precision .
integration is carried out on the 
 gpu in a manner analogous to that on the cpu since it is intrinsically 
 parallel given that each atom can be handled independently .
the integration 
 is done on the gpu rather than on the cpu to negate the need for expensive 
 copy operations of the coordinates back and forth between cpu and 
 gpu . as mentioned above ,
the berendsen and langevin thermostats are applied 
 in the same way during the coordinate update step , while the anderson 
 thermostat is implemented as a separate kernel that randomizes the 
 velocities at regular intervals .
the reason for including the berendsen 
 and langevin thermostats within the integration kernel is that they 
 operate on every time step , while the anderson thermostat is typically 
 called every few hundred time steps . in the case of the langevin thermostat , 
 a total of 3n random numbers
are required on each 
 step , which are extracted from a pool of gpu generated random numbers , 
 stored in global memory , which is refilled as needed using the parallel 
 random number generator implemented within the curand library from 
 the cuda toolkit .
the same pool is used for the anderson thermostat 
 as needed . in the interest of performance 
 and conserving limited gpu resources , the energy is only calculated 
 alongside the forces when needed .
calculation of the energies causes 
 a small degree of register spillage to global memory as well as additional 
 code execution , and since one only periodically needs these values , 
 typically only every 1 to 2 ps when writing to the output file , it 
 is not necessary to calculate these on every iteration .
additionally , 
 this is one aspect of the algorithm that is allowed to be nondeterministic 
 at the level of double - precision round - off error .
we allow this because 
 the energy values have no effect on the trajectory of the simulation . 

however , it was disturbing to observe this in action in the spsp mode 
 where nondeterminism was initially allowed at the level of single - precision 
 round - off error .
this was harmless as the forces remained consistent , 
 but in the interest of avoiding confusion among end - users , with regression 
 tests showing differences for repeat runs in spsp mode , all energy 
 summation was promoted to double - precision where differences in round - off 
 of the energies are seen only once in approximately 10 printed interactions .
the reasons for this were 
 to maximize portability of the code and avoid hardware - induced complexities 
 such as the fact that inter - gpu communication as implemented in cuda 
 4.0 requires all gpus to be on the same pcie controller .
support for 
 gpu to gpu copies in different nodes is also not yet sufficiently 
 mature to be exploited in a widely used software package . for these 
 reasons ,
the current parallel gpu scaling should be considered as 
 a lower bound on performance that will likely improve considerably 
 in subsequent versions of amber running on next generation hardware . unlike the cpu mpi implementation , the gpu mpi implementation is 
 fully deterministic for a given number of nodes and gpus . at present
, 
 our gpu implementation performs significant load - balancing between 
 the sms within each gpu rather than between gpus , which makes it possible 
 to produce a deterministic implementation . achieving good parallel 
 scaling on cpu clusters has always required the use of extensive load - balancing 
 due to the noise introduced from the operating system sharing the 
 cpu resources on a node .
gpus on the other hand provide significantly 
 more stable performance , and thus load balancing between gpus is not 
 as critical . in the current version of the software ( amber v11 ) , 
 the gpu parallel 
 support for gb calculations is implemented by dividing the force calculation 
 evenly and linearly across all gpus .
for m gpus running 
 a simulation consisting of n atoms , gpu i calculates all forces for atoms ( i  1 ) 
  n / m + 1 to i  n / m .
since gb is a full o(n ) calculation , this introduces 
 a small degree of redundancy at the gpu level with a worst - case outcome 
 of doing twice the overall calculation , but usually much less than 
 this as long as n  m. as 
 currently implemented , there are three stages to the calculation in 
 parallel , each in need of synchronizing force data across all gpus . 

the procedure used consists of three sequential mpi_allgather operations 
 per iteration to merge born radii , born force , and general force data . 

as the hardware evolves and gpus can reliably communicate directly 
 with each other within a node , we intend to scrap the internode mpi 
 communication and instead replace it with direct peer to peer copies 
 between gpus .
to assess both the serial 
 and parallel performance that can be 
 achieved with gpu accelerated amber , we ran a series of md simulations 
 representative of typical research scenarios using the gpu and cpu 
 implementations on a variety of hardware .
the systems used consisted 
 of partially folded trpcage ( 304 atoms ) , 
 ubiquitin ( 1231 atoms , pdb code 1ubq ) , apo - myoglobin 
 ( 2492 atoms ) , and nucleosome ( 25  095 atoms , pdb code 1kx5 ) . in all three 
 simulations , the ff99sb version of the 
 amber force field was used with a time step of 2 
 fs and bonds to hydrogen atoms constrained using the shake algorithm . 

the hawkins , cramer , and truhlar gb model ( amber input igb=1 ) was used with no cutoff applied to the nonbonded 
 interactions and a cutoff of 15  for the calculation of the 
 effective gb radii .
the output and trajectory files were written to 
 every 1000 steps ( 2 ps ) .
the 
 software base used for all simulations was amber version 11 , 
 including patches 115 for ambertools and patches 117 
 for amber , which were released on august 18 , 2011 .
the executables were built under the redhat enterprise 
 linux 5 operating system with intel compiler version 11.1.069 , the 
 intel mkl library version 10.1.1.019 , and the nvidia cuda compiler 
 version 4.0 .
mvapich2 version 1.5 was used for the parallel runs for 
 both the cpu and gpu versions of the code .
while the fortran compiler and use 
 of mkl has little impact on the gpu code performance , this compiler / library 
 combination gives the best performance for the cpu code and thus forms 
 a fair basis when assessing the capabilities of the gpu implementation . 

both the cpu and gpu simulations were run on machines equipped with 
 a supermicro x8dtg - df motherboard with dual hex - core intel x5670 processors 
 clocked at 2.93 ghz and mellanox quadruple data rate ( qdr ) infiniband 
 interconnects .
each compute node was equipped with one gpu and one 
 infiniband card , which were each connected to a pcie x16 slot .
error - correction 
 code ( ecc ) was switched off on the tesla cards , and the 64-bit nvidia 
 linux driver version 275.21 was used .
additionally , cpu simulations 
 were run on the xsede trestles supercomputer at the san diego supercomputer 
 center , whose nodes are equipped with four oct - core amd opteron 6136 
 processors clocked at 2.4 ghz and also interconnected with qdr infiniband .
the results of the simulation 
 timings for single gpu runs are summarized in table 2 . for small simulations such as trpcage with 304 atoms where 
 the gpu s arithmetic hardware
can not be fully utilized , the 
 serial gpu implementation is approximately twice as fast as a current 
 state of the art cpu node .
the real advantage of the amber gpu implementation , 
 however , becomes apparent for implicit solvent gb simulations of medium 
 to large systems with 2500 to 25  000 atoms .
for apo - myoglobin 
 ( 2492 atoms ) , the serial gpu version of the code is significantly 
 faster on a single gpu than a state of the art cpu node .
for nucleosome , 
 this performance gap increases dramatically such that a simulation 
 that takes two weeks on a single desktop machine with one gpu would 
 take nine months using the same desktop machine without gpu acceleration . for details on the hardware and 
 software stack ,
the 
 cpu code requires > 10 atoms 
 per core , and thus the trpcage simulation was run on 24 cpu cores .
the parallel performance of 
 the gpu and cpu implementations is compared in table 3 .
a throughput of up to 135.1 ns / day for apo - myoglobin and 
 up to 3.95 ns / day for nucleosome can be achieved with eight nvidia 
 m2090 gpus , which is a factor of 4 to 5 faster than the maximum throughput 
 that can be achieved on a typical supercomputer given that the cpu 
 scaling plateaus long before it reaches the performance achieved by 
 the gpu implementation .
even a single gpu is a factor of 2 to 3 faster 
 than the cpu scaling limit . for details on the hardware and 
 software stack ,
as mentioned above , the 
 amber gpu implementation supports three different precision models . 

the dpdp model is logically equivalent to the traditional dp approach 
 used on the cpu .
however , the use of dp on gpus can be more detrimental 
 to performance than its use on cpus .
therefore , the desire to achieve 
 high performance prompts for use of sp where possible .
it is critical 
 though to use sp carefully in order not to detrimentally affect the 
 accuracy of the md .
for this reason , we designed our hybrid precision 
 spdp model to achieve performance very close to that achievable with 
 sp but in a manner that does not impact accuracy . while we recommend 
 using the spdp precision model , it is instructive to compare the performance 
 of the different precision models available in the amber gpu implementation . 
 as shown in table 4 , the single precision ( spsp ) 
 model achieves the highest performance as expected but , as highlighted 
 in section 5 , at the cost of accuracy . our 
 hybrid spdp precision model , however , achieves a performance of greater 
 than 75% of the spsp precision model with accuracy comparable to that 
 of the full double precision ( dpdp ) model , which is between 4 to 7 
 times slower . for details on the hardware and 
 software stack ,
a critical and often overlooked 
 aspect of any major change to a 
 widely used scientific software package is a detailed validation of 
 any and all new approximations made . in the case of the amber gpu 
 implementation
the first is a detailed testing of the code itself to ensure that 
 it correctly simulates the existing systems and does not introduce 
 any new bugs or contain critical logic errors .
this is the key reason 
 why we implemented a full double precision ( dpdp ) model within the 
 gpu code .
this precision model matches the cpu code to machine precision 
 and thus can be used to check that the gpu code is indeed giving the 
 correct answers .
this allowed us to take the cpu - based regression 
 tests that we have used for many years to validate an amber installation 
 and run them using the gpu code and obtain answers to the limits of 
 machine precision .
more complex , however , is the second category for 
 validation which is the careful testing of our hybrid ( spdp ) precision 
 model .
this is significantly more complex since it requires a careful 
 evaluation of any subtle differences in the dynamics as well as ensemble 
 properties from converged simulations run in dpdp and spdp precision . 
 in this section ,
we attempt to extensively validate our spdp model 
 comparing a number of key observables across all three , dpdp , spdp , 
 and spsp , precision models .
the first metric 
 tested was the effect of the numerical precision model on the force 
 calculations as compared to a reference precision model , in this case , 
 the result from the cpu implementation .
we are using the starting 
 structures of the test systems described in section 4 above .
the dpdp model matches the reference forces very 
 closely with maximum deviations not exceeding 10 kcal/(mol  ) and rms deviations not exceeding 10 kcal/(mol  ) , even for systems as large as nucleosome .
these 
 deviations are entirely due to the different order of execution of 
 the floating point operations in the cpu and gpu implementation .
the 
 dpdp gpu implementation of pmemd will thus generate trajectories of 
 precision equivalent to the cpu implementation .
the forces obtained 
 from the spsp model , however , show very large deviations from the 
 reference values on the order of up to 10 kcal/(mol 
  ) for the largest system studied .
such large deviations introduce 
 significant numerical noise into the simulations , and it is reasonable 
 to expect that this may render the results of md simulations questionable . 
 calculating the force contributions in sp and accumulating them in 
 dp , however , leads to an improvement of more than 1 order of magnitude 
 in all cases for the spdp model as compared to the spsp model . in 
 the following section , we will indeed show that the forces obtained 
 from the spsp model are not precise enough to conserve energy in biomolecular 
 md simulations and can lead to instabilities in long time scale md 
 simulations , while the spdp model is sufficient for this purpose .
one of the most 
 important gauges for judging the precision of md software is its numerical 
 stability , which is reflected in its ability to conserve the constants 
 of motion , in particular the energy . to this end , we have performed 
 constant energy md simulations of the first three test systems described 
 in section 4 above .
we collected data for 
 100 ns ( trpcage ) , 50 ns ( ubiquitin ) , and 20 ns ( apo - myoglobin ) after 
 an initial equilibration for 1 ns at 300 k using langevin dynamics . 

in addition to simulations using a time step of 2 fs with bonds 
 to hydrogen atoms constrained using the shake algorithm with a relative 
 geometrical tolerance of 10 , we have run simulations 
 using time steps of 0.5 and 1.0 fs without constraints .
the 
 energy drifts along the trajectories are summarized in table 6 , while figure 3 shows a 
 plot of the total energy for the trajectories of trpcage and ubiquitin 
 for the different precision models at a time step of 0.5 fs .
the corresponding 
 plot for apo - myoglobin is very similar to that for ubiquitin and can 
 be found in the supporting information along 
 with plots for the larger time steps .
the plots underline that it 
 is important to validate the precision models for trajectories that 
 are long enough to uncover numerical instabilities .
while the spsp 
 model seems to perform reasonably well for a trajectory of 1 ns length , 
 in particular for smaller systems like trpcage , it becomes apparent 
 that the errors introduced lead to unacceptably large energy drifts 
 in the long term .
apart from the magnitude in the energy drift , the 
 plots for the other time steps show a similar behavior ( see supporting information ) .
table 6 shows that the spdp model is able to conserve energy to the 
 same degree as the full double precision model ( dpdp ) and the reference 
 cpu implementation .
the spsp model , on the other hand , shows energy 
 drifts that are 1 to 2 orders of magnitude larger and , on the time 
 scales investigated , leads to completely unphysical results . while 
 the cpu implementation and the gpu implementations with the full double 
 precision model ( dpdp ) and mixed precision model ( spdp ) achieve excellent 
 energy conservation on the order of 10 kt / ns / dof , 
 where dof means degree of freedom , for simulations with a time step 
 of 0.5 fs , the energy drift for the spsp model remains as large as 
 with a time step of 1.0 fs .
in addition , the energy drift for the 
 spsp model is not constant and in general increases significantly 
 along the trajectory as the temperature increases , which eventually 
 leads to unstable simulations .
total energy ( kcal / mol ) along constant 
 energy trajectories using 
 a time step of 0.5 fs without constraints . shown
are results for trpcage 
 ( left ) and ubiquitin ( right ) for different precision models .
the shake algorithm to constrain 
 bond lengths to hydrogen atoms was used for a time step of 2.0 fs ; 
 no constraints were used for smaller time steps .
a dash indicates 
 that the system heated up extremely during the simulation to the point 
 that it is meaningless to report an energy drift .
an asterisk indicates 
 that the energy drift increases dramatically for longer time scales . by way of comparison to 
 other gpu implementations of md ,
the closest published constant energy 
 gb simulation is that of the  repressor ( 1254 atoms ) using 
 openmm .
simulations of this system at 
 a time step of 1.0 fs over 1 ns of simulation time without constraints 
 using the gpu implementation of openmm gave energy drifts of 0.0054 

kt / ns / dof ( ati hardware ) ( table 
 3 of friedrichs et al . ) compared with 
 our values of 0.000006
kt / ns / dof ( dpdp ) , 0.000030 kt / ns / dof ( spdp ) , 
 and 0.000305
kt / ns / dof ( spsp ) for the similarly sized ubiquitin system 
 ( 1231 atoms ) .
it is also interesting to compare the energy conservation 
 achieved with the amber gpu implementation using the typical simulation 
 settings that one would use in a production md simulation with those 
 published for the openmm implementation .
most amber simulations would 
 use shake constraints and a 2 fs time step with the spdp precision 
 model , which for ubiquitin gives an absolute energy drift of 0.000165 
 kt / ns / dof .
this compares extremely favorably with the drift of 0.0060 
 ( nvidia hardware ) and 0.0558 ( ati hardware ) for openmm .
the energy 
 conservation of the amber gpu implementation looks even better when 
 one considers that the openmm constraint numbers were obtained with 
 a smaller time step of 1 fs . on the molecular 
 level , protein dynamics are at the heart of biological function .
md simulations that are able to accurately describe 
 protein motions aid in understanding the multitude of functions that 
 proteins carry out and can be used to interpret and predict experimental 
 results that are related to its structural and dynamic properties , 
 for example , nmr spectroscopic parameters such as spin relaxation or residual dipolar couplings . in order to 
 scrutinize the reliability of our implementation for protein dynamics , 
 we present results of md simulations of ubiquitin with the different 
 accuracy models of our gpu implementation and compare these to trajectories 
 obtained with the cpu implementation which serves as a reference . 

we focus on root - mean - square deviations ( rmsds ) and root - mean - square 
 fluctuations ( rmsfs ) of the c backbone carbon atoms 
 with respect to the crystal structure ( pdb code 1ubq(53,54 ) ) .
the highly flexible end tail of ubiquitin ( residues 71 to 76 ) 
 was excluded from our analysis .
results for the radius of gyration rg of ubiquitin can be found in the supporting information .
although the use 
 of a thermostat will to some extent cover up numerical noise introduced 
 by numerical inaccuracies in the implementation , we are analyzing 
 constant temperature simulations here since typical biomolecular amber 
 simulations are generally run with some form of temperature control . 
 in order to obtain statistically meaningful results ,
50 independent 
 md trajectories each of 100 ns length were generated at 300 k both 
 for the cpu implementation and for each of the precision models of 
 the gpu implementation .
the ff99sb force 
 field was used for all simulations with a time step of 2 fs and bonds 
 to hydrogen atoms constrained using the shake algorithm with a relative 
 geometrical tolerance of 10 .
( amber input igb=5 ) 
 was used with no cutoff applied to the nonbonded interactions and 
 a cutoff of 15  for the calculation of the effective gb radii . 

the output and trajectory files were written to every 1000 steps ( 2 
 ps ) .
temperature control during production runs was achieved with 
 the berendsen weak coupling algorithm using a time constant of t = 
 10.0 ps for the heat bath coupling .
the initial coordinates and velocities 
 that form the starting point of the 50 trajectories were generated 
 using the cpu implementation as follows .
after energy minimization 
 for 2000 steps followed by heating and an initial equilibration for 
 1 ns at 300 k using langevin dynamics with a collision frequency of 
  = 1.0 ps , snapshots were extracted at 
 time intervals of 4 ns from a constant temperature md simulation at 
 300 k. each of these snapshots was then assigned random velocities 
 corresponding to a temperature of 10 k followed by heating and equilibration 
 to 300 k using langevin dynamics for 50 ps and removal of any center 
 of mass motion that may have been introduced . using the cpu generated 
 restart files in all cases
guarantees that any numerical differences 
 observed between the various implementations can be traced back to 
 the numerical precision of the implementation and is not an artifact 
 of different initial conditions .
the plot of the rmsd values vs time ( figure 4 ) clearly shows the stability of the md simulations using 
 the cpu 
 implementation and the gpu implementation with the full double precision 
 model ( dpdp ) and the mixed single / double precision model ( spdp ) .
the 
 value of rg also remains constant throughout 
 all simulations of ubiquitin using the cpu and gpu dpdp and spdp versions 
 of the code ( see supporting information ) , indicating that the relative compactness of the protein is maintained . 

the rg value stays around 11 to 12  , 
 typical for proteins of this size , assuming a spherical conformation , and are in good agreement with the rg reported for other native - state simulations of ubiquitin .
the rmsf values for each residue for the 50 
 native - state simulations are shown in figure 5 .
the data show an excellent agreement among the cpu implementation 
 and the gpu implementation using both the dpdp and the spdp precision 
 model with the majority of the residues remaining within 2 to 3  
 of their positions in the experimental structure and somewhat larger 
 fluctuations around residue numbers 33 , 48 , and 63 . on the time scale 

investigated , no major structural change is taking place during the 
 simulation , confirming that the protein remains within a native - state 
 ensemble .
root - mean - square deviations ( rmsds ) of the c backbone 
 carbon atoms of ubiquitin ( excluding the flexible tail , residues 7176 ) 
 with respect to the crystal structure for 50 independent trajectories 
 as obtained with the cpu implementation and the gpu implementation 
 of pmemd using different precision models .
root - mean - square fluctuations ( rmsfs ) of the c backbone carbon atoms of ubiquitin residues 7176 with respect 
 to the crystal structure for 50 independent trajectories of 100 ns 
 length as obtained with the cpu implementation and the gpu implementation 
 of pmemd using different precision models .
the situation is very different for the spsp precision 
 model of 
 the gpu implementation .
while approximately half of the trajectories 
 show a similar behavior and stability as the trajectories obtained 
 with the other precision models and the reference cpu implementation 
 of pmemd , the remaining trajectories experience a very large increase 
 in the rmsd values ( figure 4 ) and a corresponding 
 increase in the rmsf values ( figure 5 ) .
the 
 reason for this behavior is that the numerical noise introduced in 
 the simulations due to rounding errors can lead to an accumulation 
 of errors such that ubiquitin starts to unfold .
this is accompanied 
 by an increase in rg ( see supporting information ) and is also reflected 
 in the large rmsf values of the residues close to the termini of the 
 protein .
since it can not be predicted whether a simulation using the 
 spsp precision model will be stable or not , no judgment can be made 
 whether the results obtained are meaningful .
the gpu implementation 
 of pmemd should thus only be used with the dpdp and the spdp precision 
 models . on the basis of the results shown here , we recommend the use 
 of the mixed - precison spdp model since the numerical results obtained 
 are comparable with those obtained with the full double precision 
 dpdp model but at significantly reduced computational cost .
we presented 
 a framework for the acceleration of classical molecular 
 dynamics simulations on graphical processing units within the pmemd 
 program of the software package amber . in this contribution , we described 
 the efficient and highly precise implementation of implicit solvent 
 generalized born molecular dynamics simulations realizing performance 
 that exceeds that achievable on any current conventional supercomputer . 

the implementation supports serial runs on a single nvidia gpu or 
 parallel runs across multiple gpus using mpi .
performance for moderately 
 sized systems can exceed 100 ns per day , making microsecond or longer 
 simulations routine on deskside workstations .
for large systems , the 
 performance differential between our gpu code and the cpu implementation 
 increases significantly . a gb simulation of the nucleosome running 
 on just one nvidia m2090 gpu completes in two weeks what a current 
 high - end workstation would take nine months to run using cpus .
this 
 massive performance increase for minimal investment in hardware will 
 transform the traditional molecular dynamics workflow , providing a 
 platform for sustained innovation .
the implementation is transparent 
 to the end user and does not 
 require any specialized knowledge beyond that required to run the 
 traditional cpu code .
our gpu implementation supports the majority 
 of the features provided by the cpu implementation including the amber 
 and charmm families of biomolecular force fields .
we have developed a 
 hybrid precision model for floating point arithmetic termed spdp that 
 provides performance close to full single - precision calculations but 
 retains energy drifts in constant energy md simulations and structural 
 properties of proteins for long time scale md simulations equivalent 
 to what is achievable with full double precision .
energy conservation 
 with our hybrid spdp model is approximately 1 to 2 orders of magnitude 
 better than what has been published for other similar gpu implementations 
 of molecular dynamics .
simulation size is limited by the memory available 
 on gpu cards , which at the time of writing is at maximum 6 gb corresponding 
 to approximately 54  000 atoms .
the current cost benefits 
 of gpus are enticing , and this is driving 
 both soft- and hardware development at a rapid pace . in a few short 
 years
, gpu - based md codes have evolved from proof - of - concept prototypes 
 to the production level amber gpu implementation we described here . 
 with gpus becoming ubiquitous in workstations and also as accelerators 
 in high - performance computing ( hpc ) platforms , the impact of our implementation 
 on the field of molecular dynamics is broad and transformative .
a 
 gpu accelerated implementation of the particle mesh ewald ( pme ) algorithm 
 for explicit solvent calculations is also available in amber 11 as 
 described in another publication .
this , 
 along with future support for advanced features including accelerated 
 md ( amd ) , umbrella sampling , replica exchange md ( remd ) , thermodynamic 
 integration ( ti ) , polarizable force fields , and constant ph simulations 
 means the use of gpus is a viable alternative to traditional supercomputers 
 for classical molecular dynamics simulations .
 we present an implementation of generalized born implicit 
 solvent 
 all - atom classical molecular dynamics ( md ) within the amber program 
 package that runs entirely on cuda enabled nvidia graphics processing 
 units ( gpus ) . 
 we discuss the algorithms that are used to exploit the 
 processing power of the gpus and show the performance that can be 
 achieved in comparison to simulations on conventional cpu clusters . 
 
 the implementation supports three different precision models in which 
 the contributions to the forces are calculated in single precision 
 floating point arithmetic but accumulated in double precision ( spdp ) , 
 or everything is computed in single precision ( spsp ) or double precision 
 ( dpdp ) . 
 in addition to performance , we have focused on understanding 
 the implications of the different precision models on the outcome 
 of implicit solvent md simulations . 
 we show results for a range of 
 tests including the accuracy of single point force evaluations and 
 energy conservation as well as structural properties pertainining 
 to protein dynamics . 
 the numerical noise due to rounding errors within 
 the spsp precision model is sufficiently large to lead to an accumulation 
 of errors which can result in unphysical trajectories for long time 
 scale simulations . 
 we recommend the use of the mixed - precision spdp 
 model since the numerical results obtained are comparable with those 
 of the full double precision dpdp model and the reference double precision 
 cpu implementation but at significantly reduced computational cost . 
 
 our implementation provides performance for gb simulations on a single 
 desktop that is on par with , and in some cases exceeds , that of traditional 
 supercomputers . 
Introduction
GPU Programming Complexities
Overview of the AMBER Implicit Solvent GPU Implementation
Performance
Validation
Conclusions and outlook