{'fraction_non_alphanumeric': 0.03685558401068923, 'fraction_numerical': 0.014118695022825966, 'mean_word_length': 4.672098016925603, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Langevin dynamics simulations of biomolecules on graphics processors

Due to the very long timescales involved (us-s), theoretical modeling of fundamental biological processes including folding, misfolding, and mechanical unraveling of biomolecules, under physiologically relevant conditions, is challenging even for distributed computing systems. Graphics Processing Units (GPUs) are emerging as an alternative programming platform to the more traditional CPUs as they provide high raw computational power that can be utilized in a wide range of scientific applications. Using a coarse-grained Self Organized Polymer (SOP) model, we have developed and tested the GPU-based implementation of Langevin simulations for proteins (SOP-GPU program). Simultaneous calculation of forces for all particles is implemented using either the particle based or the interacting pair based parallelization, which leads to a ~30-fold acceleration compared to an optimized CPU version of the program. We assess the computational performance of an end-to-end application of the SOP-GPU program, where all steps of the algorithm are running on the GPU, by profiling the associated simulation time and memory usage for a number of small proteins, long protein fibers, and large-size protein assemblies. The SOP-GPU package can now be used in the theoretical exploration of the mechanical properties of large-size protein systems to generate the force-extension and force-indentation profiles under the experimental conditions of force application, and to relate the results of single-molecule experiments in vitro and in silico.

is implemented using either the particle based or the interacting pair based parallelization, which leads to a ∼30-fold acceleration compared to an optimized CPU version of the program. We assess the computational performance of an end-to-end application of the SOP-GPU program, where all steps of the algorithm are running on the GPU, by profiling the associated simulation time and memory usage for a number of small proteins, long protein fibers, and large-size protein assemblies. The SOP-GPU package can now be used in the theoretical exploration of the mechanical properties of large-size protein systems to generate the force-extension and forceindentation profiles under the experimental conditions of force application, and to relate the results of single-molecule experiments in vitro and in silico.

I. INTRODUCTION
Mechanical functions of protein fibers such as fibronectin, fibrin fibers, microtubules, and actin filaments, are important in cytoskeletal support and cell motility [1][2][3], in cell adhesion and the formation of the extracellular matrix [4][5][6][7], and in blood clotting [8][9][10]. Physical properties of viral capsids of plant and animal viruses [11][12][13], retroviruses [14], and bacteriophages [15,16], and the transitions between their stable and unstable states determine the life cycle of many viruses, including virus maturation, and infection of cells [17]. Understanding the microscopic origin of the unique viscoelastic properties of protein fibers and the crossover from an elastic to a plastic behavior in viral capsids, as well as the control of their mechanical response to an applied mechanical force constitute major areas of research in biochemistry and biophysics.
Yet, due to the high complexity of these systems (∼10 3 −10 5 residues) and to their large size (∼ 50 − 200nm), these experiments yield results that are nearly impossible to interpret without first having some a priori information about their energy landscape [9].
Standard packages for all-atom Molecular Dynamics (MD) simulations, such as CHARMM [24], NAMD [25], and Gromacs [26] among others, are being used to access the submolecular behaviour of biomolecules. However, because all-atomic modeling is currently limited to a 10−50nm length scale and 0.1−10µs duration [27][28][29], these methods allow for the theoretical exploration of equilibrium properties of biomolecules, and reaching the biologically important ms−s timescale becomes virtually impossible even for a small system. More importantly, to fully explore the free energy landscape underlying a biological process of interest, one needs to generate a statistically representative set of trajectories. One possibility is to carry out MD simulations on manycore computer clusters, but it requires tremendous computational resources and long CPU times. For example, it takes 800, 000 CPU hours to obtain 20 short 1ns MD trajectories for the southern bean mosaic virus, which contains as many as 4.5 million atoms, on an SGI Altix 4700 cluster [30]. These limitations exclude computations as an investigative tool in the study of a range of biological problems, such as the large deformations of protein fibers, the formation of biomolecular complexes and aggregates, and the mechanical failure of viral capsids, for which experimental data are already available, thereby rendering the direct comparison of the results of experiments in vitro and in silico impossible.
Although graphics processors have been originally designed for computationally intensive graphics rendering, they have evolved over the last few years into highly parallel, multithreaded computing devices. Recent technological advances in the throughput-oriented hardware architecture of GPUs with extremely high peak arithmetic performance, which employs IEEE floating point arithmetic, have unleashed tremendous computational power that can now be utilized in general purpose scientific applications. Unlike mainstream processor architecture, GPUs devote the majority of their logic units to performing actual calculations, rather than to cache memory and flow control. Massive multithreading, fast context switching, and high memory bandwidth have enabled GPUs to tolerate latencies associated with memory calls, and to run many computational cores simultaneously (Fig. 1). Programming tools for modern GPUs include several platforms such as ATI Stream Computing [31], NVIDIA Compute Unified Device Architecture (CUDA) [32,33], and Open Computing Language (OpenCL) [34]. CUDA, a parallel computing environment (a dialect of the C and C++ programming languages), provides a high level software platform that allows a programer to define kernels that are executed in parallel by independent computational threads. The cost of one GeForce GTX 295 graphics card (from NVIDIA) with 2 GPUs (2×240 processors) is just ∼$500 (i.e. ∼$1 per processor), which makes GPUs a cheap desktop-based alternative to an expensive computer cluster.
Because GPUs differ from CPUs in several fundamental ways, CPU-based methods for molecular simulations of biomolecules cannot be easily translated nor simply adapted to a GPU. In these methods, particle-particle interactions are described by the same empirical potential energy function (force field), and the dynamics of the system in question is obtained by solving numerically the same equations of motion for all particles. Hence, there is a direct correspondence between the SIMD (Single Instruction Multiple Data) architecture of a GPU at the hardware level and the numerical routines (software) used to follow the molecular dynamics. It is then possible to execute "single instruction", i.e. calculation of the potential energy or evaluation of forces, or generation of random forces, or integration of the equations of motion, on "multiple data" sets (for all particles) at the same time over many iterations using many Arithmetic Logic Units (ALUs) running in parallel. This makes MD simulations a natural candidate to implementation on a GPU, but not all algorithms are amenable to this architecture. For an algorithm to execute efficiently on the GPU, it must be recast into a data-parallel form with independent threads running the same instruction stream on different data. There exist preliminary versions of standard packages for MD simulations of proteins implemented on a GPU, such as NAMD [28,35,36], Gromacs [37], and other applications [38][39][40]. Yet, there are no GPU-based implementations of Langevin dynamics simulations that we are aware of. In this paper, we develop and test such an implementation. Because the topology and the overall structure (geometry), rather than the atomic details, govern the force-driven molecular transitions in protein systems, we employ a coarse-grained description of proteins [41][42][43] using a Self Organized Polymer (SOP) model [44,45].
Implementing the Langevin dynamics algorithm on a GPU requires detailed understanding of the device architecture. In the next Section, we review key architectural traits of modern GPUs.
The methodology for the GPU-based implementation of Langevin simulations is presented in Section 3, where we describe the particle based and the interacting pair based parallelization approaches to force computation as well as the numerical routines for generating pseudorandom numbers, constructing the Verlet lists, and integrating forward Langevin equations to the next time step. For purposes of presentation and to focus on the most essential computational aspects, we have simplified the prosentation of the formalism as much as possible. A comparative analysis of the results of CPU-and GPU-based simulations of the mechanical unfolding for a test system represented by the all-β-strand W W domain is performed in Section 4, where we also assess the accuracy of the numerical integration. We discuss the results of the GPU-based computations in terms of the simulation time, memory usage, and computational speedup (CPU time versus GPU time) for a range of proteins including small proteins, such as the W W domain, the Ig27 domain from human titin, the C2A domain from human synaptotagmin (Syt1), the γC chain and the double-D fragment of human fibrinogen (F b), single-chain models of fibrin fibers (F b monomer and dimer), and large-size protein assembly (viral capsid HK97). The main results are summarized in Section 5.

II. NVIDIA GPU ARCHITECTURE
The generally used CPUs have most of their logical elements dedicated to cache and flow control to employ complex computational logic and to provide computational cores with fast memory access (Fig. 1). This makes a CPU capable of performing computations following a sequential workflow. On a GPU, a large number of logical elements are devoted to actual computations, and cache and flow control are reduced to a small unit. These features enable GPUs to achieve high arithmetic unit density and to perform the same computational procedure(s) simultaneously for all particles in a system by using many independent threads of execution to run parallel calculations on different data sets ( Fig. 1) [32,33]. For example, on a CPU, vector addition in M dimensions is performed in a loop, where the components of the resulting vector are computed one after another. On a GPU, this procedure can be performed using M independent threads, each of which computes just one component. Hence, on a GPU, a vector sum is computed at a cost of one addition, whereas on a CPU this cost is multiplied by M.
Yet, unlike computer clusters, where each core is capable of following its own computational protocol, contemporary compute-oriented GPUs are based on the SIMD architecture, which mandates that an identical instruction stream be fed to the large array of processing units [32].
Hence, to achieve top performance on the GPU, one has to organize a computational task into a data-parallel form with many independent threads performing the same operation(s), but on different data. In addition, the task should be compute-intensive so that, most of the time, the GPU performs actual computations rather than reading and writing data [33].
A GPU contains several multiprocessors, each having its own flow control and cache unit.
Modern GPUs from NVIDIA (Tesla C1060, GeForce GTX 285) host as many as 30 multiprocessors, each with 8 ALUs (Fig. 1) [32], which is a total of 240 computational units per chip. The GPU sits on an expansion board with its own (global or device) memory -512MB on the more conventional GeForce cards and up to 4GB on Tesla cards -separate from the DRAM accessible by the CPU. Although the GPU memory systems provide bandwidths of over 100GB/s, GPU computational cores have a low amount of cache memory, which necessitates the use of a coalesced global memory read. Because all the 240 ALUs access the global memory simultaneously, the number of memory invocations (per ALU) should be minimal to optimize the GPU performance. For example, on graphics cards from NVIDIA, the size of the cache memory is ∼20KB per multiprocessor, which is sufficient to store the local data. However, when an ALU attempts to access a large window of addresses stored in an on-board memory, the cache memory may run low causing computational cores to spend more time waiting for data (latency), rather than performing computations. With a sufficient number of threads multiplexed on each ALU, latency is effectively hidden and the GPU device achieves a high performance level without the need for a cache. For example, 30, 000 threads of execution may run concurrently on the Tesla C1060 graphics card with a peak arithmetic rate of about 900 single precision gigaflops. Of the devices currently available, graphics processors produced by NVIDIA offer the most advanced as well as the most user-friendly environment. For this reason, we used the graphics card GeForce GTX 295 (NVIDIA), which has two GPUs, each with 240 1.24GHz ALUs (30 multiprocessors) and 768MB of memory. Each multiprocessor has 16KB of cache, of which 8KB is the constant memory cache and 8KB is the texture memory cache (these numbers vary with GPU make and model).
Unlike programming on a CPU, data driven programming on a GPU entails finding a way to partition the computational problem into many identical subroutines, to improve the memory access pattern, and then to optimize the numerical procedures themselves. At a software level, a programmer can define scalar program fragments, called kernels, that are executed concur-rently in a thread block in parallel warps of 32 threads. A thread, an elementary unit of the computational workflow identified by a thread index, uses data to execute instructions that are listed in a source code, and saves the result(s) to a unique location to avoid memory conflicts.
In a single multiprocessor, 32−512 threads can be combined into a thread block. Each thread block executes a stream of instructions within a single multiprocessor and many independent blocks may be run concurrently across the entire GPU device in a grid. The software model of the GPU also includes global, shared, and local memory [32]. Threads use their own local memory to store temporary variables, and threads in a block can communicate through the fast shared memory in the same multiprocessor, but threads in different blocks cannot communicate without using the global memory, which is particularly costly in terms of clock cycles required.
Although threads in the same block can be synchronized, there is no means of synchronization of threads across several blocks. A particular global memory region (software), which corresponds to the GPU on-card memory (hardware), can be accessed by all the threads running on a GPU and by a CPU, but global memory calls should be cached or combined due to latency. Minimizing the number of memory calls can be achieved by using the shared memory, data exchange among threads in one thread block, and an optimal alignment of the computational threads with the data arrays in the global memory (coalescent memory access) [32,33]. Caching the global memory can be realized by employing texture references or constant memory.

III. LANGEVIN DYNAMICS SIMULATIONS ON A GPU
In this Section, we describe the particle based and the interacting pair based methods for the parallel computation of potentials and forces due to binary particle-particle interactions.
We also outline the numerical procedures involved in the generation of the Verlet lists and the random forces, and in the numerical integration of the Langevin equations of motion. We designed the algorithm to process as many computational fragments simultaneously on a GPU as possible in order to move the GPU into full production efficiency and to minimize the GPU/CPU communication. The algorithm decomposition, along with the workload division between the GPU and the CPU, is diagrammed in Fig. 2, where we also summarize the computational workflow on a CPU and on a GPU, including the operation on the data files for the molecular topology, the particle energies and coordinates, and the data flow between the CPU DRAM and the GPU global memory (host-GPU data transfers).
In Langevin simulations of biomolecules, molecular forces are usually described by two-body (pair) potentials, such as the harmonic potential, the FENE potential [46], the Lennard-Jones potential, etc., which differ in their mathematical form. Hence, the same generic algorithm can be employed to compute these potentials. In the pseudocode listings below, the i-th residue for a protein of N residues has a unique index i ∈ [0, N − 1]; for an array of particle coordinates r, r[i] denotes coordinates of the i-th particle. We also use the following notations: global memory reads are represented by ⇐, saving data to the GPU global memory is shown by ⇒, ⇇ denotes cached texture memory reads, and ← represents local, shared, or constant memory invocation, or assignments of variables. Consider the following computational procedure for evaluating forces on a CPU: Algorithm 1: Calculation of pairwise forces (CPU implementation).
1. P ← total number of pairs in a system for a given potential .i {index of the first particle in a particle pair} 4. j ← pairs[p].j {index of the second particle in a particle pair} There are two main optimization strategies that allow one to avoid this situation. In the first approach, all the forces for one particle are computed in one thread, which requires running N threads to obtain the force values for all particles. We refer to this procedure as the particle based parallelization approach. The use of this approach results in the same force, acting on the i-th and j-th particles, but computed twice in the i-th and j-th threads [40]. Following a different strategy, which we refer to as the interacting pair based parallelization approach, force calculations are performed for all pairs in parallel using P independent theads, and 2P force values are saved to different locations in the GPU global memory. We pursued in detail both optimization strategies, which exploit the data-parallel aspects of GPU based computing.

A. The particle based parallelization approach
In this approach, N independent threads run on a GPU concurrently, each computing all the pair potentials for each particle and summing all force values (except for the random force) to obtain the total force. Although the force acting on the i-th and j-th particles is computed twice, the number of global memory calls is reduced by a factor of 2N, and the time spent on recalculating the same potential is compensated by the time saved by not waiting to write the force data to and to read the data from the GPU global memory.
Algorithm 2: Calculation of pairwise forces using particle based parallelization.
df ← force( r, par) 11. f i ← f i ± df 12. end for 13. Output: f i The array P p keeping track of the number of pairs for all residues and the matrix P airsMap of all particle pairs are pre-generated on a CPU and fetched to the GPU global memory. P p is an N-dimensional vector of integers 0, 1, 2, . . .. Each element of this vector corresponds to a single particle, and the i-th integer value is the number of particles interacting with the i-th particle. P airsMap is the N×M matrix, where M is a maximum number from array P p . The i-th row of the matrix P airsMap corresponds to the i-th particle and contains the indices of all particles interacting with the i-th particle and the constant parameters for the potential energy function. Data in the P airsMap can be easily re-arranged for coalescent memory reads.
Coordinates of the second particle for each pair are accessed at random, which allows one to take advantage of the texture reference in global memory reads. Force values are computed in parallel by N threads as follows. First, the i-th thread reads the coordinates r i using a cached texture reference (line 3) and the number of particle pairs P p [i] (line 4). By cycling through all the pairs, formed by the i residue (lines 5−12), the thread reads the index j and the coordinates r j of the j-th particle and the constant parameters (par, lines 6−8). These are used to compute the force increment df (line 10), which is added to the resulting force f (line 11) for the i-th particle. The parameters par depend on the potential used. For example, for a covalent bond, described by a harmonic potential V H =K sp ij (r ij − r 0 ij ) 2 /2, par contains the equilibrium distance r 0 ij and the spring constant K sp ij . Verlet lists: In molecular simulations, the information about the covalent bonds and the native interactions (array P p and matrix P airsMap), obtained from the PDB structure of a protein, does not change. However, the information about nonbonded long-range interactions, describing the gradual attraction and hardcore repulsion between pairs of atoms, needs to be updated from time to time. This is the most computationally demanding component of the algorithm, since the complexity of the calculation is O(N 2 ). A common approach is to take advantage of the fact that long-range interactions vanish over some distance. This allows one to use pair lists that include pairs of particles that are closer than the cutoff distance (Verlet lists) [47]. In the particle based parallelization approach, the array P p and the matrix P airsMap have to be regenerated on a GPU in order to accelerate the computation of the potential energy using Verlet lists. This can be done by rearranging the pseudocode for particle based parallelization (Algorithm 2): Algorithm 3: Calculation of forces using particle based parallelization and Verlet lists.
Cycling over p p includes all possible residue pairs to identify pairs that are within the cutoff distance. The interparticle distances are computed in line 8. The number of particles p i that are within the cutoff distance is counted (line 11), and a newly found pair is added to the matrix P airsMap at the (i, p i )-position, i.e. copied from the matrix P ossibleP airsMap (the map of all pairs) to the matrix P airsMap (the map of pairs). Once the cycle is completed, p i is saved to the array P p [i], which stores the numbers of the residue pairs in the Verlet list.

B. The interacting pair based parallelization approach
To avoid computing the two-body potentials on a GPU twice, one can design a different computational algorithm, where each thread calculates a single pair potential for two coupled residues. Then, forces acting on the interacting particles in opposite directions are computed only once, the force values obtained are saved to different locations in the GPU global memory, and all the forces exerted on each particle are summed up to obtain the total force. This approach requires additional memory calls and a gathering subroutine for the force summation, but it enables one to accelerate simulations when the number of residues N is of the same order of magnitude as the number of ALUs, and/or when the computation of pair potentials is expensive. In the following pseudocode for the force calculation, P (number of threads) is equal to the number of interacting pairs for just one potential energy term, and each thread computes forces for one pair of residues: Algorithm 4: Calculation of pairwise forces using interacting pair based parallelization. 3. par ⇐ P airsP arameters[p] {parameters for one pair of residues p} 4. i ← pair.i {the i-th particle in the pair} 5. j ← pair.j {the j-th particle in the pair} 6. shif t i ← pair.shif t i {position in array of forces for the i-th particle} 7. shif t j ← pair.shif t j {position in array of forces for the j-th particle} {saving force for the j-th particle} Each thread identifies one pair potential using the thread index p, and reads the information about the potential from vectors pairs and P airsP arameters about the constant parameters (par), the particle identity (i and j) and particle coordinates ( r i and r j ), and the global memory addresses for saving the force values (shif t i and shif t j ). In the array F , the force values are saved to the position defined by the particle index i or j (i-th or j-th row) and by the parameters shif t i or shif t j (columns). In the array pairs, each output position i and shif t i for the i-th particle has to be unique so that each force value computed is saved to a different address in the GPU global memory. This allows one to avoid memory conflicts, but requires an additional gathering kernel for summing all the forces for a given particle obtained in the array F (in Algorithm 4): Algorithm 5: Gathering kernel for the force computation.
1. i ← GPU thread index {same as particle index} 2. f i ← 0 {resulting force due to one potential energy term} 3. P i ⇐ P p [i] {number of particle pairs for the i-th particle} Verlet lists: In the interacting pair based parallelization approach, generating a Verlet list surmounts to forming the vector pairs of all residue pairs for one potential energy term. On a GPU, constructing this vector is a formidable task, since the exact position in the list, to which information about the next residue pair should be saved, is not known. One possibility is to use the atomicAdd(...) routine from the CUDA Software Development Kit [32], which allows one to add integers in the GPU global memory without running into memory conflicts even when many threads attempt to access the same memory address at the same time. However, when many threads run in parallel, identifying new pairs and saving them one after another may result in a Verlet list that is not ordered according to the particle index. This may result, in turn, in an inefficient utilization of the cache memory. To obtain an ordered Verlet list, it has to be sorted or updated on a CPU. It is more efficient to compute interparticle distances on a GPU, copy them to the CPU DRAM, and then generate a new list.   [50,51], which produces pseudorandom numbers by generating a sequence of bits from a linear recurrence modulo 2, and forming the resulting number by taking a block of successive bits. In the space of binary vectors, the i-th element of a vector is constructed using the following transformation: y i =a 1 y i−1 +a 2 y i−2 +. . .a k y i−k , where a i are constant coeffitients. Given the initial values, y 0 , y 1 , . . . y i−1 , the i-th random integer is obtained as x i = L j=1 y is+j−1 2 −j , where s is a positive integer and L=32 is the machine word size. When a k =a q =a 0 =1, where 0<2q< k, and a i =0 for 0<s≤k−q<k≤L, the algorithm can be simplified to a series of binary operations [51]. Statistical qualities of pseudorandom numbers produced using the taus88 algorithm are high [49] when taus88 is combined with the LCG. When the periods of several components of a generator are co-prime numbers, the period of a combined generator is the product of the periods of all the components. A similar approach is used in the KISS algorithm, which combines the LCG, Tausworthe and two multiple-with-carry generators [52]. We used constant parameters which resulted in periods of 2 31 − 1, 2 30 − 1, and 2 28 − 1 for three Tausworthe generators and a period of 2 32 for the LCG generator; the period of the combined generator is ∼2 121 >2×10 36 . In the pseudocode below, superscripts h and d denote the host (CPU) and the device (GPU) memory, respectively; a section of the code, executed on the GPU, is in the same listing as the code for the CPU. In the CUDA program, the corresponding code for the GPU is organized into a separate kernel: b ← (((y 1 ≪ c 11 ) XOR y 1 ) ≫ c 21 ) 7.  Output mult × ( XOR y 1 XOR y 2 XOR y 3 XOR y 4 ) 14. end for{generating the next random number} 15

End of the GPU code
In this listing, b is a temporary unsigned integer variable, y 1 , y 2 , y 3 , and y 4 are unsigned integer random seeds for three Tausworthe generators (lines 6 to 11) and one LCG (line 12), and XOR is a binary operation of exclusive disjunction; mult=2.3283064365387×10 −10 is a multiplier that converts a resulting integer into a floating point number (between 0 and 1); c 11 =13, c 21 =19,

D. The numerical integration kernel
On a GPU, the Langevin equations of motion can be solved simultaneously for all N particles in N threads working in parallel. When the particle based parallelization is utilized, the subroutines for the force computation can be incorporated into the integration kernel. This allows a programmer to use coordinate variables, stored locally in the GPU global memory, that are read only once at the beginning of the computational procedure and are passed to the next subroutine. Since all the interactions are more or less local, texture cache can be used as well to access the coordinates in the GPU global memory. When the interacting pair based parallelization is employed, the force computations can be performed in a separate kernel and the summation of all the forces (gathering kernel, Algorithm 5) can be done inside the integration kernel. Using one kernel for the force computation, the force summation and the numerical integration minimizes the number of kernel invocations on the CPU, thus, saving time for context switching on the GPU.
force exerted on the i-th particle due to several pair potentials} 4. g i ← (g x , g y , g z ) {3D vector with 3 normally distributed random numbers} {saving coordinates to global memory at the end of each step} Each thread computes the displacement vector for just one particle. Once particle coordinates are retrieved from the GPU global memory via texture reference (line 2), they are used in the computational procedures that follow (the total force is computed in line 3). The corresponding forces f i,v , where the index v is running over different potential energy terms, are computed using ether the particle based or the interacting pair based parallelization approach. When the former approach is used, the entire computational procedure from the force computation to the random force generation, and to the numerical integration can be organized into a single kernel.
In the latter case, an additional gathering kernel is needed to compute the total force (line 3).
Continuing, particles are shifted to their new positions (line 5), which are saved to the GPU global memory (line 6). Since these coordinates are used at the next time step t n+1 , they have to be moved from the time layer r[i](t n+1 ) to the time layer r[i](t n ) at the end of each iteration.

A. The SOP model
We employed the methodology for the GPU-based realization of Langevin dynamics to develop a CUDA program for biomolecular simulations fully implemented on a GPU. To describe the molecular force field, we adapted the Self Organized Polymer (SOP) model (SOP-GPU program) [44]. Previous studies have shown that the SOP model describes well the mechanical properties of proteins, including the Green Fluorescent Protein [54] and the tubulin dimer [55].
In the SOP model, each residue is described using a single interaction center (C α -atom). The potential energy function of a protein conformation V , specified in terms of the coordinates {r}=r 1 , r 2 , . . . , r N , is given by (1) In Eq. (1), the finite extensible nonlinear elastic (FENE) potential V F EN E describes the backbone chain connectivity. The distance between two next-neighbor residues i and i+1, is r i,i+1 , while r 0 i,i+1 is its value in the native (PDB) structure, and R 0 =2Å is the tolerance in the change of a covalent bond (first term in Eq. (1)). We used the Lennard-Jones potential (V AT T N B ) to account for the non-covalent interactions that stabilize the native state (second term in Eq. (1)). We assumed that, if the noncovalently linked residues i and j (|i − j|>2) are within the cutoff distance R C =8Å, then ∆ ij =1, and zero otherwise. We used a uniform value for ǫ n =1.5kcal/mol, which quantifies the strength of the non-bonded interactions. All the non-native interactions in the V REP N B potential are described as repulsive (third term in Eq. (1)). Additional constraint are imposed on the bond angle formed by residues i, i+1, and i+2 by including the repulsive potential with parameters ǫ r =1kcal/mol and σ i,i+2 =3.8Å, which determine the strength and the range of the repulsion. To ensure the self-avoidance of the protein chain, we set σ=3.8Å (last term in Eq. (1)).

B. Benchmark simulations
We carried out test simulations of the mechanical unfolding for the all-β-strand domain W W from the human P in1 protein (PDB code 1PIN, Table I) using the SOP-GPU program. The rationale behind choosing this protein as a test system is two-fold. First, the W W domain is of particular interest to the field of protein folding and dynamics, and several research groups have expended considerable efforts to characterize the biophysical and biochemical properties of this protein [29,56,57]. Secondly, this is the smallest known independently folding all-β domain and the all-β protein architecture is the primordial structural state that can be studied experimentally using single-molecule force spectroscopic techniques such as AFM, and laser and optical tweezers [58,59]. For these reasons, the W W domain has been extensively used in the theoretical exploration of protein folding and unfolding.
We consider the following principal sources of error: (1) precision issues arising from the differences in single precision (GPU) and double precision (CPU) IEEE floating point arithmetic, where G i is the random force, and f (r i )=−(∂V (r i )/∂r i ) is the total force due to the covalent CPU versus σ f * ≈6.58pN on the GPU due to the small sample size (Fig. 3). The magnitude of the critical force for unfolding is well within the 60pN−200pN force range observed for mostly β-strand single domain proteins [58,61].

C. Accuracy of the numerical integrators
Langevin simulations fully implemented on a GPU enable one to obtain long trajectories of protein dynamics generated over as many as 10  we used an exactly solvable model of a Brownian particle X(t) evolving in a one-dimensional harmonic potential, V (X)=K sp (X − X 0 ) 2 /2, where X 0 is the equilibrium position and K sp is the molecular spring constant [62].
Under the non-equilibrium conditions of the time-dependent force application, f ext (t)=r f t, the average particle position (the end-to-end distance), computed theoretically, is given by where τ =ξ/K sp is the characteristic timescale. In pulling simulations, the average particle position at the step n+1, X(t n+1 ) , where t n =nh, can be obtained recursively from the average position obtained at the previous n-th step, X(t n ) , using the first-order integration scheme, or the second-order integration scheme, Hence, Eq. (3), (4), and (5) can be used to assess the accuracy of the numerical integrators.
We carried out calculations of ∆X(t) at room temperature for n=10 9 iterations using X 0 =0 as initial condition, which corresponds to the initial molecular extension of ∆X(0)=0, K sp =40pN/nm, which is within the 20−50pN/nm range of values observed in the experimental force-extension curves of proteins. We used the time steps of h=1, 25 and 50ps. We set κ=10pN/nm and ν 0 =1µm/s, which translates to r f =10 −5 pN/ns. These are typical values of a cantilever spring constant and a pulling speed used in AFM experiments. We set the diffusion constant to D=k B T /ξ=1.5×10 −11 cm 2 /s, which corresponds to the force-driven ∼150nm extension of the fibrinogen molecule observed in AFM experiments over time t=0.1s [21]. The slow force-driven "diffusion of the molecular extension", ∆X(t), described by the Brownian particle model with D≈10 −11 cm 2 /s, should not be confused with the free Brownian diffusion of protein molecules in aqueous solution, for which D is in the 10 −6 −10 −8 cm 2 /s range.
We found that the average extensions ∆X(t n+1 ) sim , calculated using the first order (Eq. (4)) and the second order (Eq. (5)

D. Performance measurements
We have compared the overall performance of an end-to-end application of the SOP-GPU program with the heavily tuned CPU-based implementation of the SOP model (SOP-CPU program) in describing the Langevin dynamics of the domain W W at equilibrium (Fig. 4, Table   I). To fully occupy the GPU resources, we profiled the computational performance of the SOP-GPU program as a function of the number of independent trajectories running concurrently on a single GPU. We refer to this as the many-runs-per-GPU approach. Alternatively, we could have assessed the performance of the program by running one trajectory on a single GPU, but for a range of systems of different size N, which we refer to as the one-run-per-GPU approach.
The results obtained indicate that for a small system of 34 residues (W W domain, Table I), the use of a single GPU device allows one to accelerate simulations starting from 3 independent runs (for small systems, there is insufficient parallelism to fully load the GPU), which is also equivalent to running one trajectory on a single GPU for a system of ∼10 2 residues, such as the domains Ig27 and C2A (Table I) this translates to s≈100 (Fig. 4). This implies that the use of small thread blocks is more advantageous when M th ≤3, 000, e.g., when simulating one trajectory for a system of the size of fibrinogen monomer F b (N≈2, 000) or 30 runs for a system with N≈100 (Ig27), or ∼50 runs for a smaller system such as the W W domain (N=34). However, larger thread blocks should be used when M th >3, 000 to simulate a few trajectories for larger systems such as the fibrinogen dimer (F b) 2 (N=3, 849) or to obtain one trajectory for a very large system, e.g., the viral capsid HK97 (N=115, 140, see Table I).
To profile the associated computational time and memory demand for the SOP-GPU software we ran a series of benchmark simulations on several different protein systems (Table I) using the one-run-per-GPU approach. The execution time of an end-to-end application of the program as a function of particle count was recorded using the standard CUDA runtime profiling tool. To more closely assess the performance characteristics of the SOP-GPU code as a function of the system size N, we analyzed one simulation run for each protein system, generated over n=10 6 steps of size h=40ps and block size of B=64 threads (Fig. 5). For all test systems of less than ∼3, 000 residues the associated simulation time remains roughly the same, i.e. (1 − 3)×10 −4 seconds per step. This is not surprising since M th ≈3, 000 is the amount of threads needed to fully utilize the GPU resources. For larger systems (N>3, 000), the computational time scales roughly linearly with N (Fig. 5). Small variation in runtime versus N around the monotonic linear dependence is due to the different native topology of the test systems, i.e. the number of native and non-native contacts in the PDB structure (Table I) the developed formalism can be used in conjunction with more sophisticated biomolecular force fields to explore, e.g., protein-protein and protein-DNA interactions, and it can also be extended to include the three-body (angle) potentials to describe side chains. In addition, the numerical integration kernel can be modified to follow the Langevin dynamics in the underdamped limit, in order to describe the thermodynamics of biomolecular transitions.
The developed formalism has been mapped into a standard CUDA code for Langevin simulations of biomolecules (SOP-GPU program). Benchmark simulations have shown that for a test system of the all-β W W domain the results of simulations of the mechanical denaturation on the GPU agree well with the results obtained on the CPU (Fig. 3). In a separate work, we also compared the forced unfolding data, obtained on the CPU and on the GPU, for the human synaptotagmin Syt1 (manuscript under review) and for the human fibrinogen F b molecules (manuscript in preparation), and found that the results of the CPU-and GPU-based computations agree very well. Using an exactly solvable model of a Brownian particle evolving in a harmonic potential, we have assessed the accuracy of the numerical integration of Langevin equations of motion in describing the force-induced elongation of a protein chain. We found that using the first-order integrator is sufficient to accurately describe the force-driven elongation of a protein over many billions of iterations.
GPUs can be utilized to generate a few trajectories of Langevin dynamics for a large system of many thousands of residues (one-run-per-GPU approach) or many trajectories for a small system composed of a few hundreds of amino acids (many-runs-per-GPU approach). Using the onerun-per-GPU approach, we analyzed the relative CPU/GPU performance of the program, and found that the GPU-based realization leads to a substantial ∼25−30 computational speedup, which depends on the number of threads in a thread block (similar acceleration can also be achieved using the one-run-per-GPU approach). We profiled the SOP-GPU program in terms of the computational time and the memory usage for a range of proteins (Table I). Our results show that the simulation time on a GPU remains roughly constant for a system of N≈10 2 −10 3 residues, and scales linearly with N for a system of N>10 3 residues. The GPU on-board memory in contemporary graphics cards (GeForce 280, GeForce 295, and Tesla C1060) is sufficient to describe Langevin dynamics for large systems involving as many as ∼10 5 residues (Fig. 5).
Describing mechanism(s) and multiple pathways underlying biomolecular transitions and resolving the entire distributions of the molecular characteristics requires gathering of the statistically significant amount of data. This can be achieved on a single GPU using the many-runs-per-GPU approach for smaller systems, and on multiple GPUs employing the one-run-per-GPU approach for larger systems. The many-runs-per-GPU approach can also be utilized in conjunction with parallel tempering algorithms, including variants of the Replica Exchange Method, to resolve the phase diagrams of biomolecules. Beyond that, we note that due to rapid evolution of GPU hardware, the simulation time on the GPU will decrease significantly with the introduction of the MIMD (Multiple Instruction Multiple Data) based Fermi architecture (NVIDIA) in the near future [63]. The CUDA program, tested in this work, should be able to run with minor modifications on future hardware making use of increased processing resources, and will permit the theoretical exploration of a range of interesting problems for which the experimental data are already available. The presented formalism can also be applied to the studies of other systems, including molecular motors, colloidosomes and liposomes. Most of the GPU chip surface is devoted to computational units, and a graphics card has onboard global memory. On the GPU, ALUs are grouped into multiprocessors, each of which has its own small flow control and cache units. To execute a simulation protocol on the GPU device, the CPU and the GPU communicate via the PCI Express bus, which allows the CPU to execute computational kernels on the GPU and to access GPU global memory. CUDA kernels involves the decomposition of work into small fragments that can be mapped into thread blocks, and further decomposition into warps and into independent threads of execution.
The computational workflow for only one (i-th) thread running on the GPU and the workload division between CPU and GPU are shown in detail. The execution of the program is initiated on the CPU, which is used to prepare and store the initial and output data. The CPU device starts the launch of each computational kernel on the GPU for the calculation of forces, generation of the random forces and Verlet lists, and for the numerical integration of the Langevin equations of motion, using many independent threads running in parallel.  (Table I) obtained on a CPU and on a GPU (the color code is explained in the graphs). Panel    HK97). The information about the native state topology for these biomolecules is summarized in Table I.  k PDB code for a structural unit; the full HK97 capsid structure can be found in the Viper [64] database.
l Based on a cut-off distance of 200Å.