{'fraction_non_alphanumeric': 0.04094485000050624, 'fraction_numerical': 0.01457976854617433, 'mean_word_length': 4.442962636393696, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Protofold II: Enhanced Model and Implementation for Kinetostatic Protein Folding

A reliable prediction of 3D protein structures from sequence data remains a big challenge due to both theoretical and computational difficulties. We have previously shown that our kinetostatic compliance method (KCM) implemented into the Protofold package can overcome some of the key difficulties faced by other de novo structure prediction methods, such as the very small time steps required by the molecular dynamics (MD) approaches or the very large number of samples needed by the Monte Carlo (MC) sampling techniques. In this article, we improve the free energy formulation used in Protofold by including the typically underrated entropic effects, imparted due to differences in hydrophobicity of the chemical groups, which dominate the folding of most water-soluble proteins. In addition to the model enhancement, we revisit the numerical implementation by redesigning the algorithms and introducing efficient data structures that reduce the expected complexity from quadratic to linear. Moreover, we develop and optimize parallel implementations of the algorithms on both central and graphics processing units (CPU/GPU) achieving speed-ups up to two orders of magnitude on the GPU. Our simulations are consistent with the general behavior observed in the folding process in aqueous solvent, confirming the effectiveness of model improvements. We report on the folding process at multiple levels; namely, the formation of secondary structural elements and tertiary interactions between secondary elements or across larger domains. We also observe significant enhancements in running times that make the folding simulation tractable for large molecules.


Introduction
Proteins are large biomolecules that are responsible for a vast array of biological functions inside the cell, and appear in the form of enzymes, antibodies, motor proteins, transport proteins, etc. [3]. The function of a protein strongly depends on its 3D structure (i.e., 'conformation') which in turn can be directly determined from the linear sequence of amino acids (AAs) linked together to form the protein chain (i.e., 'configuration') [4]. 1 Therefore, the computer-aided prediction of the folded structure of a protein from the knowledge of its sequence (referred to as 'protein folding') is the key to understanding many biological processes in the cell. This knowledge is crucial toward the ultimate goal of modeling proper function or malfunction at molecular and cellular level (e.g., deadly diseases such as cancer, Alzheimer's, Parkinson's, etc.) and is central to a variety of bioengineering applications including 'protein engineering' [5,6].

Related Work
There are several different computational approaches for protein folding prediction, ranging from knowledge-based techniques to methods starting from physical principles [7].
Knowledge-Based Methods. The knowledge-based approaches predict the structure of a given protein using the information extracted from previously determined structures and known types of folds. They are generally more reliable than physics-based methods, but have limited applicability in predicting new types of folds. Examples of knowledge-based techniques are homology or comparative modeling [8][9][10] and fold recognition or threading [11][12][13]. We refer the reader to [7] for a comprehensive review of such methods.
Physics-Based Methods. On the other hand, methods that utilize models formulated empirically or obtained from physical principles are less reliable and more time-consuming, but apply to a wider range of folding simulations [7]. These methods range from de novo [14,15] to ab initio [16,17] prediction techniques. Here we briefly review some of the common ab initio approaches, namely sampling methods and MD simulations [7].
Sampling methods generate numerous samples in the conformation space, followed by an evaluation of their free energies. Different search algorithms are used to find the unchallenged global minimum of the free energy, assumed to be associated with the native structure according to the 'thermodynamic hypothesis' [4]. These search methods include simulated annealing [18,19], basin hopping [20][21][22][23], evolutionary algorithms [24][25][26] and MC simulation with biased moves [27][28][29]. A review of conformation sampling methods for protein folding can be found in [30]. Sampling methods have two major limitations, namely 1) they do not provide any information about the biological pathway; and 2) finding the global minimum is not guaranteed because of the finite number of samples.
On the other hand, MD approaches simulate the biological pathway using a model built upon physical principles. Standard MD techniques include Newtonian dynamics [31][32][33][34], Langevin dynamics [35][36][37][38] and Brownian dynamics [39][40][41][42]. A review of MD simulation methods for protein folding is provided in [43]. In order to keep the numerical algorithms stable, very small time steps (in the order of femtoseconds) along the simulation trajectory are required, which does not support folding simulation of typical proteins that span milliseconds except for small molecules [43].
Kinetostatic Compliance Method. The KCM was introduced in [44,45] to overcome some of the key challenges in the aforementioned approaches. In this method, implemented in the software package Protofold [46][47][48], the protein chain is modeled as a kinematic linkage which complies under the kinetostatic effect of the force-field obtained form intramolecular interactions between the atoms. The key contributions of KCM were 1. modeling the constrained kinematics of the protein chain with significantly fewer degrees of freedom (DOF) than, for example, those of the 'beads and springs' model used in many MD methods; and 2. converging faster to the minimum energy state by using kinetostatic (i.e., 1st-order) variations rather than dynamic (i.e., 2nd-order) response.
In KCM, each rotatable joint, used to model the constrained motion of the chemical bonds, changes by an amount proportional to the effective torque on that joint. It was shown that KCM is a faster and more stable alternative to the traditional dynamic simulation techniques [48]. The Protofold platform has since provided a kinematic testbed for subsequent research activities. Examples are predicting hydrogen bond connectivity sub-graphs [49], its application to the design of stable peptide nano-particles [50], the analysis of protein mobility (using the 'pebble game' algorithm) [51], the development of mechanical models for secondary structural elements [52], and nano-machanism synthesis from a 'link soup' of pre-specified structural elements [53,54].
In the earlier stages of the development, the energetics were limited to intramolecular interactions in the gas-phase of the protein-e.g., Coulombic and van der Waals forces exchanged between atoms of the protein itself, ignoring the interactions with solvent molecules. However, an important class of biologically significant proteins are water-soluble, whose folding process is predominantly driven by the interactions with the solvent, particularly the so-called 'hydrophobic effect' 2 which was missing from Protofold I [48].
Computing Solvation Effects. From a computational perspective, the solvation effects can be modeled in a number of different ways, broadly classified into 'explicit' and 'implicit' techniques.
The explicit methods use all-atom force field models such as SPC, SPC/E, TIP3P, TIP5P [55,56], or coarse-grained (CG) models [57,58] which are less structured representations of the solvent obtained by mapping two or more particles onto a single interaction site [55]. A prohibitive computational cost is associated with the large number of solvent molecules required to model a bulk solution.
Alternatively, approximate schemes that include the solvent effects implicitly can provide useful quantitative estimates, yet remain computationally inexpensive [59]. The implicit models estimate the contribution of each solventexposed atom to the solvation free energy using empirical formulae, most commonly expressed as a linear function of the solvent accessible surface area (SASA) [60]. An exact computation of SASA requires obtaining the surface area of the envelope of overlapping spheres [61], which is computationally expensive. Alternatively, approximate formulations have been developed to efficiently predict the expected (i.e., probabilistic average) SASA based on simplifying assumptions on the distribution of the coordinates of atoms (or groups of atoms) in the 3D space. For instance, CHARMM [62] uses the probabilistic approach from [63], which estimates the SASA as a function of the distances between pairs of atoms or residues. A similar model with similar parametrization [64] was used in GROMOS [65], a recent improvement of which was given in [66]. AMBER [67,68] uses a fast linear combination of pairwise overlaps (LCPO) algorithm [69], which improves the method in [63] by adding more terms to the approximation. Although being widely popular in well-known MD packages, these methods rely on simplifying assumptions that compromise accuracy. For example, the method in [63] assumes a uniform random spatial distribution of atoms or residues, which introduces bias into the simulation results.
The inclusion of the solvation free energy computed using an adequately accurate evaluation of the SASA [1,2] results in a more realistic energy-and force-model for simulating the natural behavior of water-soluble proteins.

Outline & Contributions
In Section 2, we introduce an improved free energy model making use of the linear implicit model given in [60] to compute the solvation free energy-and force-field from a knowledge of the SASA and its gradient for individual atoms at a given 3D conformation. We develop a simple offset surface enumeration technique that can approximate the SASA and its gradient up to any desired accuracy. Our method is significantly more accurate than the probabilistic methods such as those given in [63,69] yet noteably faster than the exact method given in [61], while the trade-off between speed and accuracy can be adjusted by a proper choice of the enumeration (i.e., surface sampling) density.
A second major contribution of this work is to develop significantly more efficient algorithms and data structures in Section 3 to speed up the computations, and to implement them into Protofold II: 1. We use a 3D hash table data structure based on a uniform spatial grid that supports fast queries to identify pairs of proximal atoms. This helps speeding up the computations by eliminating negligibly small interactions associated with pairs of atoms that are farther than a so-called 'cut-off distance'.
2. We use a tree-based data structure to span the protein chain efficiently for the purpose of characterizing interaction types between pairs of atoms based on their distances along the topological structure of the chain.
3. We develop sequential and parallel surface enumeration algorithms to compute the SASA and its gradient for individual atoms needed for the solvation energy and force computations, respectively, up to desired accuracy. 4. We employ prefix sums [70] to compute the joint torques on the kinematic linkage of the protein chain, given the resultant forces on the individual atoms.
As a result, the numerical complexity for each KCM cycle, including the computation of resultant forces on the atoms and the links (except those resulting from solvation effects) and their conversion to joint torques, is reduced from O(n 2 ) in Protofold I [48] to O(n) in Protofold II, where n is the number of atoms in the protein molecule. 3 The SASA evaluations for solvation force computations in our model turns out to be the bottleneck to the entire simulation-up to several orders of magnitude slower than the electrostatic and van der Waals force computations. Fortunately, the surface enumeration algorithm lends itself well to high-throughput data parallelism. In Section 4 we first present the CPU-parallel implementation using OpenMP, leading to moderate speed-up factors (up to one order of magnitude). To leverage the massive processing power offered by the single-instruction multiple-thread (SIMT) architecture of the modern graphics hardware-onto which our data-parallel SASA enumeration algorithm maps perfectly-we present a GPU-parallel implementation and its optimization. The implementation takes advantage of the device memory hierarchy and hiding memory access latencies, in turn leading to larger speed-ups (up to two orders of magnitude).

Formulation
Section 2.1 starts with an overview of the underlying kinematic principles of the KCM simulation first introduced in [44][45][46][47][48]. The protein chain is modeled as an open kinematic linkage with reduced DOF in terms of dihedral and rotamer angles, which complies under the effect of interatomic and solvation forces. Next, the energy-and force-field formulation used in Protofold II is described in Section 2.2, with special emphasis on the newly introduced solvation effects. Lastly, the KCM optimization process is presented in Section 2.3.

Kinematic Model
Proteins are long polymeric chains made of AAs, which exist in only 20 different types (except for few rare exceptions), joined together as a linear polypeptide chain [3], structural details of which are summarized in Appendix A. Here we restrict ourselves to the kinematic representation of the chain's conformation within the scope of KCM.
Linkage Parameterization. Figure 1 schematically illustrates the repetitive sequence of −N−C α −C− atoms, 4 called the 'backbone' or the 'main chain', with 'side chains' resembling branches that extend out of it. As explained in Appendix A, the backbone conformation can be specified to an adequate accuracy by two sets of dihedral angles; namely, where m is the number of AA residues along the chain. The conformation of each side chain, on the other hand, can be specified by up to 4 extra angles −180 • ≤ χ i,k < +180 • for 1 ≤ k ≤ l i where 0 ≤ l i ≤ 4 is the number of side chain links of AA i , and the subscript k corresponds to the bonds numbered in the obvious order along the side chain C and N atoms.
To set a reference for the angle measurements, the zeroreference position description (ZRPD) method [71] is used. The zero-position (ZP) for the protein chain is defined as the conformation of the serial linkage in which all peptide planes are coplanar (i.e., φ 0 i = ψ 0 i = −180 • ) and side chain dihedrals are set to default low energy values identified as 'rotamers' [45]. To unify the notations, all angular variables are denoted by where 0 • ≤ θ j,k < 360 • . The shifts in (1) and (2) by the intercept values of φ 0 i = ψ 0 i = −180 • and in (3) by the favorable rotamer angles χ 0 i,k ensure θ j,k = 0 • at the ZP conformation. A similar indexing scheme is used to identify the unit vectors along the rotation axes of revolute joints associated with these angles denoted by u j,k (1 ≤ j ≤ 2m, 0 ≤ k ≤ 4), i.e., • u 2i−1,0 (1 ≤ i ≤ m) is the unit vector along the bond between N of AA i and C α of AA i ; is the unit vector along the bond between C α of AA i to the C of AA i ; and are the unit vectors along the successive side chain C and N atoms.
Thus the kinematics of the linkage-which abstracts the protein conformation-can be completely specified in terms of the rigid body rotation transformations obtained from these rotation angles and rotation axes. The spatial orientation of the rigid peptide planes can be described conveniently with a pair of base vectors whose linear combination spans the peptide plane. The so-called 'body vectors' are denoted by b j,k (1 ≤ j ≤ 2m, 0 ≤ k ≤ 4), i.e., is the base vector that connects the C α of AA i to the N of AA i+1 ; and are the base vectors along the successive side chain C and N atoms.
The first two sets of vectors listed above are called the 'main chain body vectors'. Every vector in the peptide plane that describes the relative positions of any two atoms can be obtained as a linear combination of these base vectors as C 1 b 2i,0 + C 2 b 2i+1,0 . The coefficients C 1 and C 2 , referred to as 'peptide plane constants', are invariant with respect to the rotations in the chain, thus can be precomputed prior to the KCM simulation. Different pairs of coefficients are used for vectors describing the relative positions of different pairs of atoms. Based on experimental evidence, it is a reasonable assumption that these coefficients are the same across all AAs [72], and the average values are given in Table 1. 5 In addition to the main chain body vectors, the 'side chain body vectors' (the third group listed above) are defined for the relative positions of the C and N atoms along the side chains. We refer the reader to [72] for more details about the molecular model of the peptide unit.
For a protein chain with m AA residues, the number of links can be obtained as noting that l i ≤ 4. The term 2m counts two rigid links per each AA's backbone-one for −CO−NH− and one for −C α −  [72]. in the peptide unit-in order to have each rigid link connected to the next with a single revolute joint along either N−C α or C α −C, as depicted in Fig. 1. The second term accounts for the number of additional side chain links. As a result, the total DOF of the kinematic linkage is equal to the number of links. Table 2 gives a complete description of dihedral angles, unit vectors, and body vectors for the entire protein chain.
Kinematic Equations. The instantaneous conformation of the protein chain is related to the ZP conformation with a set of rigid body transformations. Given a particular conformation in terms of θ j,k (1 ≤ j ≤ 2m, 0 ≤ k ≤ 4), the unit vectors and body vectors are transformed as follows: where the superscript 0 indicates the reference ZP conformation. [M j,k ] is the 3 × 3 matrix representation of the rigid body transformation M j,k ∈ SO(3) that maps the ZP unit and body vectors u 0 j,k and b 0 j,k to their transformed orientations u j,k and b j,k , respectively. These vectors are expressed using 3 × 1 column matrices. The transformation matrix for the main chain vectors (k = 0) can be calculated as a product of successive rotations around individual joints in the main chain: while the transformation matrix for the side chain vectors (k ≥ 1) is defined as a product of rotations around joints in the main chain, and those around the side chain: where [R r,s ] is the 3 × 3 matrix representation of the joint rotation transformation R r,s ∈ SO(3) around the ZP unit vector u 0 r,s with an angle θ r,s (1 ≤ r ≤ i, 0 ≤ s ≤ k) [45], using the right-hand rule to choose the direction.
Once the body vectors are obtained using (6) for a given conformation, the moved atom center positions can be computed for the individual atoms. Assuming that the N atom at the amino-terminus is fixed at the origin, the coordinates of the main chain N and C α atoms are obtained as where the index j = 2i − 1 corresponds to the C α atom of residue AA i while the index j = 2i corresponds to the N atom of the residue AA i+1 for 1 ≤ i ≤ m. The coordinates for the other atoms in the peptide group, namely H, C and O, are obtained from those for C α and N, and a linear combination C 1 b 2i,0 + C 2 b 2i+1,0 of main chain body vectors using the coefficients C 1 and C 2 given in Table 1. For the side chain C and N atoms, the coordinates are similarly obtained as where k = 1, 2, 3, and 4 corresponds to the successive side chain C and N atoms in the residue AA i . The coordinates for all other side chain atoms are obtained similarly from vectors along the side chain bonds subjected to the same set of motions [44].
The atom position vectors obtained from (9) and (10) at each iteration are used in Section 2.2 to compute the energies, forces, and torques that will determine the motion for the subsequent iteration.

Force Model
The interatomic effects can be classified into covalent and noncovalent interactions. The covalent interactions need not be considered explicitly in the force-field, since they are implicitly introduced in terms of the kinematic constraints innate to the kinematic chain model adopted in Section 2.1. The noncovalent forces, which are responsible for conformational changes in the protein chain, can be derived from the free energy formulation that follows.
For a protein chain made of n atoms, we assign each atom with a unique index 1 ≤ i ≤ n, and its center coordinates r i ∈ R 3 obtained from dihedral angles using kinematic equations (9) and (10). 6 Each atom is identified by a unique tuple a i = (i, r i , R i , q i , i , γ i , · · · ) (1 ≤ i ≤ n), containing its index, position, radius, charge, well depth parameter, solvation parameter, and other atomic constants, to be introduced shortly. The set of all the atoms in the molecule is denoted by A = {a 1 , a 2 , · · · , a n }. The aggregated free energy of all the atoms in A can be decomposed into the following terms: where G elec (A) is the electrostatic energy, including intramolecular charge interactions, hydrogen bonding effects, and the induced polarization in the solvent when the molecule is dissolved. G vdw (A) is the sum of intramolecular van der Waals energies, also called 'steric effects', resulted from induced dipoles in the molecule. The sum of the first two terms has been accounted for in Protofold I [46][47][48] using the AM-BER force-field model [67]. G cav (A) is the nonpolar solvation free energy, the free energy change resulting from transfering a molecule from vacuum to solvent, i.e., the entropic change due to the formation of the cavity occupied by the instantaneous 3D shape of the protein [73]. Experimental results have shown that many water-soluble protein folding reactions are predominantly driven by a favorable reduction in ∆G cav (A) [3], hence we incorporated this term into the improved energy formulation for Protofold II.
Electrostatic Interactions. The charge interactions are formulated using the modified form of Coulomb's law [67]: where d i,j = r i − r j 2 is the interatomic center distance, q i and q j are the electrostatic charges, and r i and r j are the position vectors of the pair of atoms a i , a j ∈ A, respectively. ε i,j = κ i,j ε 0 is the 'dielectric constant' and is generally larger than vacuum permittivity ε 0 ≈ 8.854 × 10 −12 Farads (i.e., κ i,j > 1). Thus (12) can be used to calculate the interactions between charges in the solvent, if a continuum dielectric model is used [3]. The dielectric constant reflects the ability of the environment to attenuate electrostatic interactions, e.g., κ i,j ∼ 80 for aqueous solvent and κ i,j ∼ 2-4 for the interior of the protein [3], where the larger value for the former is due to the induced polarization of water molecules. It is worthwhile noting that because of the nonuniformity of the dielectric medium, the most accurate computation of the electrostatic energy requires solving Poisson-Boltzman (PB) differential equations [74]. However, solving PB for every cycle of the KCM simulation is computationally expensive, and approximate alternatives such as generalized Born (GB) model can be used [75,76]. A simple distance-dependent dielectric constant is used here (following the convention in [48]) to mimic the polarization effect, with closer interactions weighted more heavily [67]. The resultant Coulombic force F elec i = −∇ ri G elec applied on the atom a i by other atoms is then obtained as where e i,j = (r i − r j )/d i,j is the unit vector along the line of centers of the pair of atoms a i , a j ∈ A. Since F elec electrostatic interactions between atoms that are farther than a so-called cut-off distance d elec cut := 9.0Å are usually deemed negligible in the literature [3]. 7 Therefore (13) is approximated as follows to reduce the number of pairwise computations between all the atoms: where A elec } is referred to as the neighborhood of atom a i associated with the electrostatic force-field, and consists of all the atoms whose distance to a i are bounded by the cut-off distance d elec cut .
Van der Waals Interactions. The van der Waals interactions are formulated using the empirical Lennard-Jones 6-12 potential function formula [67]: where d i,j = r i − r j 2 is the interatomic center distance, D i,j = R i + R j is the 'van der Waals distance' in which R i , R j are the van der Waals radii of the atoms a i , a j ∈ A, respectively. i,j = √ i j is the 'depth of potential well' for the particular pair of atoms. It is worthwhile noting that the van der Waals effects have the same origin as the electrostatic forces, and reflect the induced dipoles due to transient fluctuations in electron clouds of the interacting atoms [3]. The resultant van der Waals force F vdw i = −∇ ri G vdW on the atom a i by other atoms is then obtained as where e i,j = (r i − r j )/d i,j is the unit vector along the line of centers of the pair of atoms a i , a j ∈ A. The van der Waals forces have a much smaller effect radius and are significant only when the atoms approach each other very closely. The repulsive component becomes very large when the two atoms penetrate into each other, an effect widely known as the 'steric clash'. The attractive component known as the 'London dispersion' force, on the other hand, is dominant when the atoms are farther than the van der Waals distance D i,j [3]. These interactions decay much faster than Coulombic forces, hence a smaller cut-off distance of d vdw cut := 5.0Å is sufficient [3] resulting in the following approximation: where } is referred to as the neighborhood of the atom a i associated with the van der Waals force-field, and consists of all the atoms whose distance to a i are bounded by the cut-off distance d vdw cut . Interaction Classification. The interactions discussed so far are between the pairs of atoms that are not covalently bonded, thus (14) and (17) have to be modified to eliminate the terms corresponding to the pairs of bonded atoms (i.e., '1-2 interactions'). Furthermore, it is a common convention in molecular mechanics to modify the electrostatic and van der Waals interactions between the pairs of atoms bonded to a common atom, i.e., atoms that are 2 bonds apart along the chain (i.e., '1-3 interactions'), as well as the atoms that are 3 bonds apart along the chain (i.e., '1-4 interactions') [77]. Consequently, the empirical forms of (14) and (17) are modified as where w elec i,j and w vdw i,j are the weight factors for the electrostatic and van der Waals forces, respectively, for the pair of atoms a i , a j ∈ A depending on their interaction type. The weights are set to w i,j = 0 for 1-2 interactions, and 0 ≤ w i,j ≤ 1 for 1-3 and 1-4 interactions, whose values vary across different force models [62,65,67,68]. w i,j = 1 for all other situations. In other words, the atoms that have at least 4 bonds in between them along the graph of covalent bonds are far enough to be considered unaffected by the covalent electron clouds, as originally formulated in (14) and (17).
Nonpolar Solvation Effects. The hydrophobic free energy of solvation, which reflects the entropy changes in the solvent molecules due to cavity creation, is formulated using the linear empirical formulation in [78]: where γ j is the atomic solvation parameter and A exp j is the Lee and Richards SASA for the atom a j ∈ A [81]. To ob- Table 3: Atomic solvation parameters adopted from [78] for different adjustments by [79,80].
Atom type Adjustment in [79] Adjustment in [80] tain the atomic SASA at a given snapshot, a probe radius of R H2O = 1.2 − 1.4Å is used to offset the van der Waals surfaces of the individual atoms as illustrated in Fig. 2. These offset spheres are overlapped to obtain the contributions of different atoms to the total SASA. The atomic solvation parameter γ j is an estimate of the free energy per unit area required to transfer the atom from vacuum to water, and depends on the atom type [78]. Table 3 shows the values of this parameter for different atom types (namely, C, uncharged O or N, S, O − , and N + ) obtained in [78] based on the experimental results in [82] adjusted by [79,80]. The hydrophobic interaction forces F cav i = −∇ ri G cav on the atom a i by other atoms is then obtained as where ∇ ri A exp j is the gradient of the exposed area of the atom a j due to an infinitesimal displacement of a i . It is important to note that, unlike the force formulae presented earlier for the electrostatic and van der Waals effects in (13) and (16), the summation in (21) for the solvation free energy gradient iterates over all a j ∈ A including a i itself.
Considering the case when i = j, one realizes that displacing the atom a i in any direction has the same effect on the geometry of the protein surface as displacing all the other atoms in the opposite direction. Therefore Substituting for this term in (21) leads to the following symmetric form, whose range of summation excludes a i itself, similar to (13) and (16): We show in Section 4 that (23) is computationally preferable over (21). To further simplify (23), note that for a pair of atoms a i and a j the infinitesimal displacement of one does not affect the overlapped solvent exposed area of the other if their offset spheres (i.e., the spheres with radii R off H2O . Therefore, the number of terms that contribute a nonzero value to the summation of (23) is significantly reduced: where A cav } is referred to as the neighborhood of the atom a i associated with the nonpolar solvent effects. For practical reasons that will be explained in Section 3.2, we use a larger neighborhood redefined as A cav cut } using the more conservative (but constant across all pairs of atoms) cut-off distance of d cav cut : 0Å is typically safe. Note that unlike the case with (13) and (16), eliminating pairwise interactions with d i,j > 8.0Å from (24) does not impart an approximation error.

Kinetostatic Simulation
We use the KCM (presented in [44][45][46][47][48] for protein folding) to explicitly integrate the conformational changes of the linkage model under the kinetostatic effect of the force-field computed in Section 2.2.
Link Forces and Torques. For a protein chain with a total of l = O(m) links, where m is the number of AA residues, the resultant force and torque applied to the j th link (1 ≤ j ≤ l) are computed as where r i is the absolute center position vector of the atom a i ∈ A obtained from (9) and (10) in Section 2.1 (with different index notation), and L j ⊂ A is the subset of atoms that belong to the j th link along the chain. Note that the origin of the absolute coordinate system (arbitrarily picked the same as the N-terminus anchor of the chain) is selected as the torque center for all links.
Equivalent Joint Torques. Since the revolute joints are assumed to be frictionless, the action of the link forces and torques can be replaced by an equivalent set of torques acting on the joints [48], as shown in Fig. 3. For a given rigid link, one can trace a unique serial chain of h successive links (1 ≤ h ≤ l) starting from the N-terminus and ending at the link under consideration, which is equivalent to a path along the graph of the open linkage. The contribution of the force F link h and the torque T link h applied to the end-effector link (i.e., the h th link along the serial chain) to the joint torque at the k th joint along the chain preceding the end-effector, Figure 3: The forces on each link are converted into an equivalent set of joint torques on the preceding joints in the chain.
, can be computed using the conventional manipulator Jacobian matrix [J] [48]: where [η h ] = [η 1,h , η 2,h , · · · , η h,h ] T represents an h × 1 array of joint torques that the end-effector force and torque will induce on the different joints preceding the end-effector along the serial chain.
[J] T is the transpose of the 6 × h Jacobian for a given configuration of the chain [48]. The k th column of the Jacobian associated with the k th revolute joint is given by where u k is the unit vector along the k th joint and (p h − p k ) is a vector connecting a point p k along the k th joint's axis to the point p h where the end-effector force F h applies-i.e., the atom positions obtained from (9) and (10). This process is repeated for all links to calculate the contribution of each link on the joints preceding that link in the serial chain. The total torque for each joint is obtained from a summation of these contributions [48]: where the indexing h = k, k + 1, · · · , l of the links is ordered along the main chain from amino-terminus to carboxylterminus, and can branch along the side chains. The range of the summation in (29) implies that each joint torque τ k (1 ≤ k ≤ l) depends on the forces F link h and torques T link h on the succeeding links (k ≤ h ≤ l).
Kinetostatic Simulation. Making use of the assumption in [48] that the inertia forces at the atomic scale have negligible effects on the dynamics of folding compared to interatomic forces, the successive kinetostatic fold compliance relates the joint torques to the changes in the dihedral angles as follows: where ∆θ j and τ j are the angular change and the joint torque of the j th revolute joint (1 ≤ j ≤ l), respectively. |τ max | is the maximum joint torque throughout the entire chain, used to normalize the torques to the interval [0, 1], and the coefficient κ is chosen small enough to avoid large changes in the angles, and to achieve numerical stability. One can notice that the conformational change computed in (30) is analogous to taking a step along the steepest-descent direction of the free energy gradient in the conformation landscape.
The computed changes in the joint angles are applied to the kinematic chain, and the entire process is repeated on the updated conformation until a convergence criteria is reached, as described in more detail is Section 4.
It is worthwhile noting that once the chain conformation (i.e., optimization variables) is modeled as in Section 2.1 and the energy-field (i.e., objective function) is formulated as in Section 2.2, the search for local or global minima of the free energy in (11) can be undertaken using a variety of classical (e.g., conjugate-gradients) and stochastic (e.g., genetic algorithm) optimization methods. Since the focus of this article is mainly on force-field modeling and computing, we skip a detailed treatment of the search phase.

Algorithms
This section presents efficient algorithms and data structures to speed up the force field computation during kinetostatic iterations. To leverage the proximity information between the atoms, we use a 3D hash table data structure based on a uniform spatial grid in Section 3.2. To classify the interatomic interaction types based on chain topology to compute the electrostatic and van der Waals effects, we use a treebased data structure in Section 3.3 that replaces the adjacency matrix used in Protofold I [46][47][48]. To compute the solvation effects, we develop an approximate (yet adequately accurate) surface enumeration technique in 3.4, efficient CPUand GPU-parallel implementations of which will be detailed in Section 4. Finally, we compute joint torques by aggregating contributions of different links (traversed along different paths in the linkage graph) on the joints along the chain, using the well-known 'prefix computation' in Section 3.5 which can also be implemented efficiently in parallel [70]. We show that the computational complexity of all steps is decreased from O(n 2 ) in Protofold I [46][47][48] to expected O(n) in Protofold II for a protein chain with a total of n atoms.

Rigid Transformations
At every snapshot t ≥ 0 of the KCM, the protein conformation is described by a set of dihedral angles θ t j,k defined in (1) through (3).
• At the first iteration (t = 0), all angles are initialized as θ 0 j,k = 0 (ZP conformation). • At subsequent iterations (t ≥ 1), for 1 ≤ j ≤ 2m (where m is the number of AA residues) and 0 ≤ k ≤ l i (where 0 ≤ l i ≤ 4 is the number of side chain links of the residue AA i ), the angles are obtained as θ t j,k = θ t−1 j,k + ∆θ t−1 j,k , where the increment ∆θ t−1 j,k is computed using (30) from the previous step's configuration and joint torques.
Once the dihedral angles are known, the transformation matrices [M t j,k ] are obtained from (7) and (8) using sequential matrix multiplication traversing the linkage tree from the anchored amino-terminus to the open carboxyl-terminus. Next, the unit vectors u t j,k and the body vectors b t j,k are computed from (5) and (6). Since the number of links is clearly less than the number of atoms, these computations take O(n) time. The Cartesian coordinates of the individual atom center positions r i ∈ R 3 (1 ≤ i ≤ n) are obtained from the body vectors using (9) and (10), which also takes O(n). Hereon, we assume that both dihedral angles and atom center positions are known for the purpose of computing the next snapshot's energies, forces, and torques.

Proximity Queries
The brute-force approach for obtaining the proximity information at each snapshot is to check center distances against the cut-off distance for all possible pairs of atoms, which takes O(n 2 ) time. Using this method, the approximate truncated formulae given in (12), (15), and (20) would take the same asymptotic time as the exact all-pairs formulae given in (13), (16), and (23), respectively, which is O(n 2 ). Geometric hashing provides a simple solution to speed up proximity queries. Table. We use a 3D hash table data structure based on a uniform Cartesian grid, which bounds the current 3D structure of the protein and arranges the atom indices into groups based on their center positions. Each 3D grid cell is associated with a so-called 'bucket' that stores the indices of the atoms whose centers are located inside the cell in a linkedlist. The grid dimensions are set dynamically to adapt to the shape of the protein's bounding box at the current snapshot.

3D Hash
The grid cells are chosen to be cubic, i.e., with equal edge length s c along all 3 Cartesian axes. Given the min/max corner coordinates of the bounding box of the atom centers r min , r max ∈ R 3 -which can be obtained in O(n) by scanning through the n atom center coordinates-we choose s c in such a way that it results in αn grid cells/buckets, where α > 0 is an arbitrary constant. More precisely, we choose where v BB (A) is the protein bounding box volume. The dimensions of the grid bounding box are then chosen as (r max − r min )/s c s c (slightly larger than the dimensions of the protein bounding box r max − r min ), where the operator · is applied componentwise along the 3 Cartesian axes. Before we proceed with presenting the complexity analysis, we make the following assumptions: Assumption 1 Due to the extremely strong repulsive van der Waals forces, the atoms that are not covalently bonded cannot penetrate into each other, and those covalently bonded intersect over a small volume. Given any arbitrary subset of atoms A ⊆ A with R min = min ai∈A R i and R max = max ai∈A R i , let the maximum penetration volume between any pair of covalently bonded atoms a i , a j ∈ A be upper-bounded by min{v i is the volume of the atom a i , and 0 ≤ < 1 4 is a small number. Since each atom makes at most 4 covalent bonds, the unpenetrated volume for the atom a i is lower-bounded by ( Assumption 2 It is also reasonable to assume that if the protein is either in an extended conformation aligned with one of the Cartesian axes (which is the case near the ZP conformation) or in a globular conformation (which is the case for most water-soluble proteins at their folded conformation), the empty space inside the bounding box is not extremely larger than the space occupied by the protein atoms, i.e., v BB (A)/v(A) = O(1). Supported by experimentation, we assume this to be the case in the intermediate conformations as well, to simplify the analysis. However, there are possible conformations (e.g., extended along a diagonal direction in the axis-aligned bounding box) that would violate this assumption and result in slightly larger running times than predicted here, in spite of the low probability.
Letting A := A (hence |A | = n) in Assumption 1, and noting from the definition that v BB (A) = (αn)s 3 c , from Assumption 2 it follows that 3α 4π (s c /R) 3 = O(1). Therefore, if we choose the grid cell size s c ∼ 1Å then α = O(1) and the number of buckets will be αn = O(n). Table Construction. At each snapshot of the simulation, the algorithm scans through all the atoms with updated positions, and the deterministic hash function simply maps the atom center positions r i ∈ R 3 that lie inside a grid cell into the corresponding 3D array of buckets. The 3D index k ∈ Z 3 ∩ [0, +∞] 3 of the bucket to which a given atom a i ∈ A belongs is determined in O(1) time as k = (r i − r min )/s c , where the operator · is applied componentwise along the 3 Cartesian axes. Therefore, scanning through the atoms and constructing the 3D grid data structure is expected to O(n) time and to requires O(n) space.
Neighbor Queries. Once the atoms are arranged into the buckets, the algorithm iterates through the grid cells and scans through the linked-lists within the buckets. For each atom a i in a given bucket associated with the grid cell index k, the set of 'neighbor atoms' defined as can be identified rapidly for a given cut-off distance d cut associated with any of the energetic interactions explained in Section 2.2. As illustrated in Fig. 4, a spherical region of radius r c = d cut + √ 3s c is considered around the (center point of) each grid cell to look for the (center point of) 'neighbor cells', defined as the collection of cells which completely lie inside this spherical region. The cut-off distance d cut is offset by the diagonal size of the cells √ 3s c which takes into account the worst-case difference of the distance between cell centers and the distance between atom centers. This guarantees that the set of all the atoms inside this collection of covered cells (denoted as A i ) contains the set of all neighbor atoms, i.e., , sincē R, s c ∼ 1Å and d cut ∼ 10Å. As a result, it is expected to take O(1) time to scan through the atoms inside the collection of neighbor cells, and O(n) total time to traverse all pairs of atoms using the 3D hash table.
For parallel implementation purposes, we construct (and dynamically maintain) a 'neighborhood matrix' composed of an array of n O(1)−sized linked-lists (one list per atom), where the i th list (1 ≤ i ≤ n) contains the indices of the neighbor atoms A i . Constructing this data structure is expected to take  (18) and (19), respectively. It is important to note that such interactions are strictly pairwise, i.e., they depend on the relative positions of pairs of atoms a i ∈ A and a j ∈ A elec i or A vdw i , and the presence of any third atom a k ∈ A (k = i, j) does not affect the force exchanged between a i and a j . Therefore, the com-putation algorithm is straightforward: it iterates over all the atoms for 1 ≤ i ≤ n (sequentially or in parallel) and for each atom, it computes F elec i (A i ) and F vdw i (A i ), by sequentially aggregating the contributions of the hashed neighbor atoms a j ∈ A i , using (18) and (19), respectively.
Unfortunately, this is not the case for solvation force computation using (23) or (24), which requires computing the gradients of the atomic SASA with respect to the coordinates of the set of neighbor atoms. The SASA variations in one atom a i ∈ A with respect to an infinitesimal change in the position of another atom a j ∈ A cav i , can be affected by the presence of a third atom a k ∈ A cav i (k = i, j), thus cannot be obtained in a pairwise fashion. This is because the overlaps of the pairs of offset spheres are not mutually disjoint. A segment of the offset surface can be overlapped with more than one neighbor sphere simultaneously, thus displacing one of the overlapping spheres may or may not affect the SASA. We return to this subject in Section 3.4 where a surface sampling algorithm is proposed as a simple solution.

Bonds Tree/Graph
To decide the weight factors w elec i,j in (18) and w vdw i,j in (19), one needs to quickly identify the types of interactions based on the number of bonds between pairs of atoms as described in Section 2.2. An n × n look-up table was used in Protofold I [46][47][48] for all pairs of n atoms which required a preprocessing step with O(n 2 ) time and space. This can be improved by constructing a tree/graph data structure that stores the combinatorial structure of the chain, in which the vertices are atoms and edges are the covalent bonds between them. By excluding a single edge from the loops associated with the rare aromatic groups in certain side chains (e.g., imidazole in His and indole in Trp), this graph can be converted to a tree whose root is arbitrarily chosen as the N atom of the amino-terminus. The interaction types are then identified by the shortest path lengths between atoms.
The algorithm starts from the root and visits all vertices using a standard tree traversal routine. For each vertex, it stores a pointer to its parent and the index of the corresponding atom's AA residue. During the force computation in each KCM iteration, the residue indices for a pair of atoms of interest are checked. If the atoms are farther than a residue apart (i.e., if AA indices are neither identical nor consecutive), the weights in (18) and (19) are simply set to 1 (i.e., 1-4 interactions or beyond). Otherwise, the algorithm checks 1) if one atom is the parent of the other (i.e., 1-2 interaction); 2) if one atom is the grand parent or sibling of the other (i.e., 1-3 interaction); or 3) if one atom is the great-grand parent or sibling of parent of the other (1-4 interactions). This requires O(n) time and space for preprocessing and O(1) query time during the KCM iterations.

Surface Enumeration
As mentioned in Section 1.1, several attempts have been made to approximate the SASA and its derivative by a pair-wise treatment of the overlaps, including the probabilistic methods [63,64,66,69], popular in many molecular simulation software such as CHARMM [62], GROMOS [65], and AMBER [67,68]. In an early attempt to add solvation effects to Protofold II, we used these pairwise approximate formulae, which made it possible to compute the solvation forces with running times comparable to those of the electrostatic and van der Waals force computations. However, a comparison with the exact method [61] showed that when the distribution of the atoms deviates from that assumed in the the probabilistic methods, prohibitively large errors can be introduced into the resultant effects. The exact method [61] takes O(|A i |) operations for computing the SASA and its gradient for a i ∈ A by using the coordinates of all neighbor atoms A cav i ⊆ A i . This is asymptotically O(1) time per atom (since we reasoned earlier that |A i | = O(1) under Assumption 1), but in practice it is not nearly as fast as using the pairwise formulae. Alternatively, we use an approximate method that relies on an enumeration of the surface area, in which the deviations from the exact results can be controlled to a desired precision in a trade-off with computation time.
Offset Sphere Sampling. For a given atom a i ∈ A of van der Waals radius R i , an offset sphere of radius R off i = R i + R H2O concentric with the atom sphere is considered. The atom's SASA is obtained by measuring the area A exp i of the portion of the offset surface that is not overlapped by the offset sphere of any neighbor atom a j ∈ A cav i (hence exposed to the solvent). To approximate the SASA, one can generate a large but finite 'quasi-uniform' set of sample points denoted as Q i (1 ≤ i ≤ n) on the surface of the offset sphere of the atom a i ∈ A, by which we mean a sampling that allows approximating the exposed fraction of the surface by the ratio of the number of exposed sample points to the total number of sample points. In other words, if we let be the subset of the solvent-exposed sample points, i.e., the points that are outside the offset spheres of all neighbors, then lim |Qi|→∞ |Q exp using a large enough sample size |Q i | 1. 8 There are different ways to obtain a quasi-uniform deterministic sampling on a sphere with consistent incremental quality [84]. For example, one could use a triangular spherical meshing algorithm, which starts from an icosahedron approximation of the sphere and recursively creates successive triangular subdivisions projected back on the sphere. Alternatively, one could use a polar geodesic sampling algorithm, 8 An alternative approach is uniform random sampling, e.g., using the simple method in [83]. Random sampling is easier to implement in parallel since every sample point in Q i would be independent from others, and results in A exp ] is the expected ratio of the exposed sample points in probabilistic terms. However, it requires much larger sample sizes to approach the expectation and to achieve adequate accuracy in practice.
which starts from a set of orbits with uniform angular distribution and samples a number of points uniformly on each orbit proportional to the orbit's circumference. We take the latter approach whose details are presented in [2].
To improve the efficiency, one could always precompute the coordinates for a single sampling Q on a unit sphere centered at the origin, and map it into individual atoms with different offset sphere center positions r i and radii R off i using the mapping Q i = T i (Q) where T i (q) = r i + R off i q for 1 ≤ i ≤ n and q ∈ Q. This implies an equal number of sample points N = |Q| for all the atoms, selected as N = 4π(R off max ) 2 /δA where R off max = R max + R H2O is the maximum offset sphere radius, and δA is the desired characteristic area element carried by each sample point. Hence the sampling takes O(nN ) operations for all the atoms regardless of the sampling technique. For implementation purposes, we assign an arbitrary ordering to the sample points, letting Q = {q 1 , q 2 , · · · , q N } and denoting transformed sample points as q i,k = T i (q k ). (20), the total solvation free energy can be computed directly from

Energy and Force Approximations. Substituting for
To obtain the solvation force on any atom a i ∈ A, the energy must be differentiated with respect to the atom's center coordinates r i ∈ R 3 , giving rise to (23). It is very important to note that an infinitesimal displacement of the atom a i can change the SASA of a i itself, as well as that of the neighbor atoms a j ∈ A cav i ⊆ A i . However, we showed in Section 2.2 that the two effects are geometrically dependent, yielding the symmetric form in (24). From a computational perspective, (24) is preferred over (24), because 1. it eliminates the need for computing the gradient ∇ rj A exp i for the cases when i = j, hence decreasing the number of such computations from n 2 to n(n − 1); and 2. its symmetric form lends itself to a data-parallel implementation that is balanced between computation and data sharing tasks, as we show in Section 4.
This follows from the fact that for a given pair of indices i = j, the first term in the formula for F cav i is identical to the second term in the formula for F cav j (both given by (24)), cutting the number of required SASA gradient computations to n(n − 1)/2. Substituting for A exp i and A exp j in (24), the solvation forces can be approximated as where ∇ rj f exp i and ∇ ri f exp j can be approximated using the forward-difference method from finite variations of f exp i and f exp j with respect to the positions of the atoms a i and a j , respectively: where δr > 0 is the finite difference, e s (s = 1, 2, 3) are the unit vectors along the 3 Cartesian axes, and f exp i,j,s are the exposure ratio of the atom a i ∈ A after changing the position of the neighbor atom a j ∈ A i from the current value r j to a hypothetical variant r j + δre s .
Enumeration Algorithm. In order to compute the exposure ratio and its finite-difference gradient, we use a binary enumeration function B : on the offset sphere of the atom a i ∈ A is overlapped by at least one neighbor offset sphere (i.e., if ∃a j ∈ A i such that T i (q k ) − r j 2 ≤ R off j ) and B(a i , q k ) = 0 if the sample point is exposed to the solvent. The algorithm iterates over all the atoms for 1 ≤ i ≤ n and all the sample points for 1 ≤ k ≤ N (sequentially or in parallel). For each sample point, the indicator B i,k := B(a i , q k ) is initialized to 0, and each point is tested against the set of neighbors A i , scanned sequentially. As soon as one overlapping neighbor is found, b i,k is set to 1 and there is no need to test the rest of the neighbors. The exposure ratio is then computed as In the worst case, this takes |A i |N tests and N binary sums per atom where N = |Q| is the sample size, which adds to O(N ) basic operations per atom (since we reasoned earlier that |A i | = O(1) under Assumption 1), and a total of O(nN ) time for all the atoms. For every sample point, the sequential inner loop of the algorithm can be repeated 3|A i | times for computing the variations f exp i,j,1 , f exp i,j,2 , and f exp i,j,3 used in (35), after introducing the finite difference to the 3 Cartesian coordinates (one at a time) of each neighbor atom a j ∈ A i . This takes 3|A i | 2 N more tests per atom, still asymptotically O(N ) but not fast enough in practice. There is a notably more efficient way to do the latter computation by ruling out the subset of sample points that cannot possibly contribute to f exp i,j,s − f exp i (s = 1, 2, 3) in (35) during the first iteration when computing B i,k indicators. In particular, if a sample point is overlapped by more than one neighbor, displacing any neighbor does not affect its exposure state (from overlapped to exposed or vice versa), hence it does not contribute to f exp i,j,s − f exp i . To leverage this property, we expand the binary definition of the state function to C : A × Q → Z ∩ [0, ∞) such that C(a i , q k ) counts the actual number of neighbors a j ∈ A i that overlap the sample point q i,k ∈ T i (Q). Three different states for a sample point are observed in terms of the changes in C i,k := C(a i , q k ): Algorithm 1: SASA enumeration algorithm for solvation free energy and force computation.
1. 'Not overlapped' or 'exposed' (C i,k = 0). In this case, displacing any neighbor either keeps the state at C i,k = 0 or changes it to C i,k = 1, where the latter case affects the contribution to SASA. Hence the inner loop needs to be repeated for all neighbors (i.e., for 3|A i | times).
2. 'Critically overlapped' (C i,k = 1). In this case, the only neighbor whose displacement may change the sample point's state to C i,k = 0 is the one that originally overlapped it, and displacing any other neighbor either keeps the state at C i,k = 1 or changes it to C i,k = 2, both of which correspond to overlapped states that does not affect the contribution to SASA. Hence the inner loop is repeated only 3 times after displacing that critical neighbor along the 3 Cartesian axes.
3. 'Multiply overlapped' (C i,k ≥ 2). In this case, displacing any neighbor either keeps the state at C i,k ≥ 2 or changes it to C i,k = 1, both of which correspond to overlapped states that does not affect the contribution to SASA. Hence the inner loop need not be repeated at all.
Therefore, the only changes that contribute a nonzero value to f exp i,j,s − f exp i (s = 1, 2, 3) are those from exposed (C i,k = 0) to critically overlapped (C i,k = 1) and vice versa, thus a significant amount of computation time can be saved by early detection of the rest. An atom a j ∈ A i is called a 'critical neighbor' of the atom a i ∈ A with respect to a sample point q i,k ∈ Q i along a particular direction e s (s = 1, 2, 3), if a finite displacement δre s results in such a change. As a direct consequence of geometry, if a j is a critical neighbor of a i along +e s , then a i is also a critical neighbor of a j along −e s , both with respect to the same sample point. Therefore, a pair of neighbor atoms a i , a j ∈ A exchange a solvation force ±δF cav = ±4πγ i (R off i ) 2 /(N δr)e s due to their overlap at the sample point q i,k if and only if they are critical neighbors with respect to e s . The improved algorithm (based on the integer-valued C i,k ) is different from the original (based on the binary-valued B i,k ) in that the first iteration of the sequential inner loop for computing C i,k terminates after the second (rather than the first) overlap is encountered, because all C i,k ≥ 2 have equivalent implications according to the above rules. 9 During this step, the value of f exp i is initialized to 1 for each atom, and every time a sample point with C i,k = 1 or 2 is discovered, f exp i is decremented by 1/N . The next 3 repetitions of the inner loop per neighbor atom displacement depend on the aforementioned rules based on the value of C i,k . The 3 variants of the exposure ratio f exp i,j,s (s = 1, 2, 3) are initialized to f exp i for each atom with respect to displacements in all of its neighbors. Every time a sample point with C i,k = 0 or 1 is encountered, the inner loop is repeated with displaced neighbor coordinates to discover the critical neighbors, each adding ±1/N to f exp i,j,s . Significant speed-ups are achieved in terms of the average time, a rigorous analysis of which is not possible without assumptions on the spatial distribution of atoms. However, the worst case time complexity is still O(nN ) for the sequential algorithm. One could easily parallelize the algorithm at the outer loops over the atoms and sample points, while the inner loops over the neighbor atoms is best implemented sequentially. On a simple concurrent-read concurrent-write (CRCW) parallel random-access machine (PRAM) with common conflict resolution (briefly introduced in Appendix C.1), the parallel running time of O(nN/P ) can be achieved in theory using P processors (i.e., linear speed-up), which leads to O(n) if we have P = O(N ) processors at our disposal-not far from reality when using GPUs. However, there are more complications to the machine architecture in practice, as will be addressed in Section 4.
The complete process is described in pseudo-code in Algorithm 1. The instructions marked by a dagger ( † ) modify variables that belong to different atoms iterated in parallel by the outer-most loop, hence require a mutex with nuances that depend on the architecture as described in Section 4.

Prefix Sum Calls
There are multiple references in Protofold II to the generic prefix sum routine-explained in Appendix B, which can be performed using optimal sequential and parallel algorithms in linear number of steps-that emerge naturally as a consequence of the linear topology of the polypeptide backbone: 1. Computing the link transformations from the successive matrix multiplications in (7) and (8), in which the domain is Σ := SO(3) (represented by 3 × 3 rotation matrices) and the operator ⊕ is the matrix multiplication.
2. Computing the atom center coordinates from the successive vector summations in (9) and (10), in which the domain is Σ := R 3 (represented by 3 × 1 column matrices) and the operator ⊕ is the vector summation.
3. Computing the joint torques from the successive superposition of the contributions of each link on the preceding joints in the chain using (27) and (29), in which the domain is Σ := R 6 (represented by 6 × 1 column matrices) and the operator ⊕ is the inner product. The contribution of P h on the k th joint is obtained as the inner product of the two matrices η k,h = [J k ] T [P h ] arranged into the following matrix:  (27). Note that each link only affects the preceding joints in the chain, hence η k,h = 0 for all h ≤ k − 1. The total torque joints τ k (1 ≤ k ≤ l) can be obtained as a summation over the rows of the above matrix via (29). In Protofold I [46][47][48] this was accomplished by scanning through the terms along the columns in (37), which took l(l + 1)/2 = O(l 2 ) operations. In Protofold II we perform row scanning of the matrix, starting from the bottom row and moving upwards. More specifically, by factoring out the Jacobian terms [J T k ] in each row of (37)

Protofold II Architecture
Unlike Protofold I [46][47][48] that was programmed in Matlab , Protofold II is reprogrammed with a new architecture in C++. The CPU-and GPU-parallel algorithms are implemented as external modules and linked to the main application thread as dynamic link libraries (DLL), which can be integrated into other folding packages.
As depicted in Fig. 5, a typical KCM simulation in Protofold II can be summarized into the following steps: 1. Input: The user specifies a primary structure (i.e., AA sequence information) to the interface.

Preprocessing:
The program constructs the AA chain using the structural assumptions given in Section 2.1 to arrange the atoms into the consecutive peptide planes. The double-bond angles are all set to the fixed values of ω i = 0 • (cis) or −180 • (trans) and the body vectors are assigned with the values given in Table 1.

Initialization:
The main chain dihedral angles are initialized as φ 0 i = ψ 0 i = −180 • and the side chain dihedral angles are initialized to rotameric default values χ 0 i,k [45] (3) referred to as ZP initial conditions.

Forward Kinematics:
The conformation variables summarized in Table 2 are converted to the Cartesian coordinates of the individual atoms by using the sequence of rigid body transformations described in Section 2.1.

Coordinate Hashing:
Using the 3D grid data structure presented in Section 3.2, the atom coordinates are arranged into buckets for fast neighborhood queries based on the cut-off distances.
6. Force Computations: The free energy-and forcefields are computed from the atom coordinates using the equations given in Section 2.2. This is where the CPUor GPU-parallel modules are called for computing the solvation effects.

Torque Computations:
The forces on the atoms are converted to joint torques using the Jacobian transformation described in Section 2.2.

KCM Stepping:
The kinetostatic effect of the joint torques are computed using the simple steepest-descent stepping explained in Section 2.3.

9.
Termination: If the convergence criteria is met, the program terminates; otherwise it repeats the steps 4 through 8 above.
10. Output: The intermediate (every several frames) and final conformations in PDB format, the variations of the dihedral angles and free energy terms, and the performance measures (e.g., running times of different steps) are exported by the program.
These steps characterize the process of arriving from sequence configuration (i.e., primary structure) to stable 3D conformation (i.e., tertiary structure) without any additional assumptions. Although this is the ultimate goal of protein folding, it is rather ambitious to obtain results that are consistent with experimental observations except in the case of relatively short chains; e.g., folding simulation of α−helix coiling described in section 5.1. This is due to a variety of reasons ranging from the sensitivity of the folding pathway to the physical parameters (e.g., adjusted coefficients in the empirical force-field equations) to the sensitivity of the spatial structure of long chains to simplifying geometric assumptions (e.g., the exact planarity of the peptide planes).
Additional Functionalities. In order to enable addressing certain computer-aided structural studies on real proteins effectively in spite of the aforementioned difficulties, we found it imperative to include the following additional functionalities in Protofold II: • The user has the option to 1) specify only sequence data, from which the 'canonical' peptide plane geometry (i.e., assuming exact planarity ω i = 0 • (cis) or −180 • (trans) and average lengths in Table 2); or 2) import the protein structure as a PDB file and retain the peptide group geometry as-read when constructing the rigid links.
• The user has the option to limit the mobility of the linkage by fixing as many dihedral angles as desired. This enables folding studies at multiple levels and different scales. For example, it is possible to group collections of AAs (e.g., secondary elements, motifs, domains, etc.) into presumed rigid bodies and limit the DOF to deformations at the loops connecting them.
• In addition to the ZP initial conditions, the user may choose to use other initial conditions, including but not limited to completely random initial conditions or the native conformation perturbed by arbitrary (deterministic or randomized) changes to certain dihedral angles.
• When importing PDB files, the program eliminates water molecules-since their effect is implicitly incorporated by the solvation energies-but retains other heteroatoms (e.g., metal ions, co-factors, substrates, etc.) and includes them among chain atoms when computing the force-field. This is crucial since the proper folding of many proteins is dependent on these agents.
In addition to the above features, the following need to be included in future versions: • The program currently supports monomeric protein folding in its simplest topology. It is desirable to enable multimeric protein folding by maintaining multiple chains bound together (i.e., quaternary structure) and more complex topologies induced by other effects (e.g., disulfide bonds, hydrogen bonds, lipidation, etc.) • the simplistic steepest-descent search process presented in Section 2.3 has not evolved much since Protofold I [46][47][48]. Our numerical experiments suggest that better optimization algorithms such as a hybrid Monte Carlo sampling combined with steepest-descent or conjugategradients KCM 10 could be more effective in avoiding local minima and enable faster convergence to the global minimum.
Parallel Implementations. As demonstrated in Section 3.4, the solvation energy and force computations using (24) and Algorithm 1 are the most time-consuming steps of each KCM iteration, mainly due to the large number of sample points |Q| = N 1 required to enumerate the offset sphere of each atom for an adequate approximation of SASA and its gradient. To benefit from the single-instruction multiple-data (SIMD) characteristic of Algorithm 1, the variables pertaining to different atoms are assigned to different processors. The two terms on the right-hand side of (24) are computed concurrently by different processors assigned to a i , a j ∈ A and broadcasted to each other to minimize the computational work. An immediate consequence is an additional communication overhead and possible network contention due to concurrent write attempts. Such a trade-off between computation and communication intensities is a common characteristic of parallel algorithms [85], and will be considered here for code optimization.
Here we focus on the implementation of the SIMD Algorithm 1 using two parallel computing models; namely, • one that is designed for coarse-grained multiprocesser machines such as multi-core CPUs (Section 4.2); and • another that is designed for fine-grained multiprocesser machines such as many-core GPUs (Section 4.3).

CPU-Parallel Implementation
The first parallel implementation targets a coarse-grained shared-memory multiprocessor machine, i.e., one with a multi-core multi-thread CPU. Given n atoms n threads are generated, assigning one thread per atom a i (1 ≤ i ≤ n).
The neighborhood information (i.e., the list of indices of all a j ∈ A i defined in (31) for each atom a i ∈ A) is constructed and saved in the global memory shared among all processors, hence can be accessed concurrently from different threads. Each thread stores the sample point coordinates, the SASA and its gradient, and the resulting solvation energy and force components in the local memory of the processor. The thread iterates sequentially over the sample points on the offset sphere, and a counter variable that keeps track of the number of overlapped sample points is initialized within the scope of the thread. For each sample point, the coordinates are computed and tested sequentially against all neighbors to obtain Once the exposure states are obtained for the original configuration of the neighbors, the thread loops over all neighbors one more time to examine the effects of their displacement along the 3 coordinate axes one at a time. If the criteria given in Section 3.4 are met, the pair of force components ±δF cav = ±4πγ i (R off i ) 2 /(N δr) need to be added to the total solvation forces of two neighbor atoms a i and a j along the proper coordinate axis, and in opposite directions. This results in two write operations per incidence, the first of which modifies F cav i of a i , which is safely assigned to the current thread and occurs in the local memory without any concern related to communication between the threads. The second write operation, on the other hand, modifies F cav j of a j , a variable assigned to a different thread. This requires communication between the two threads, and has to be implemented using atomic write operations into the global memory to guarantee mutual exclusion. Figure 6 shows the multithreading scheme for the CPU-parallel algorithm. The algorithm is implemented using the OpenMP library. Although linear speed-up is expected in theory on an abstract CRCW PRAM, the actual speed-up is sublinear (as depicted in Section 5.2) in practice due to bus traffic, network contention, cache invalidations, and serialized operations.
CPU Optimization. The number of CPU cores is generally much smaller that the number of atoms (p n). Nevertheless, it is good practice to generate more threads than the number of cores to maximize the performance by keeping the processors saturated at all times with computational work. Accessing global memory incurs latency at the incidence of a cache miss and multithreading is a standard technique for hiding such latencies. The computation instructions are interleaved with memory access instructions, hence every time one thread is accessing the global memory the processor can switch the context to a different thread. Other optimization attempts include using local memories instead of global memories whenever possible, and avoiding multiple computations of constant parameters or variables that are used repeatedly.

GPU-Parallel Implementation
The second parallel implementation targets a fine-grained machine with a hierarchical memory architecture, i.e., one with a many-core many-thread GPU. Given n atoms, a linear grid of n blocks is generated, assigning one block per atom a i (1 ≤ i ≤ n). Each block is further divided into N = |Q| threads, assigning one thread per sample point q k ∈ Q (1 ≤ q ≤ n) on the unit offset sphere. Prior to GPU kernel execution, the neighborhood information (i.e., the list of indices of all a j ∈ A i defined in (31) for each atom a i ∈ A) is transferred from the CPU (i.e., host) memory to the GPU (i.e., device) memory. For each thread, the iteration over different neighbors is performed sequentially, similar to the CPU-parallel code presented in 4.2. The solvation energy and force components, a counter that keeps track of the number of overlapped sample points, and exposure states are initialized in the shared memory of the blocks, which require atomic operations for access safety by multiple threads, while sample point coordinates are stored in the registers that are local to each thread. Each sample point is tested sequentially against all neighbors to obtain C i,k (1 ≤ i ≤ n, 1 ≤ k ≤ |Q i |).
Once the exposure states are obtained for the original configuration of the neighbors, the thread loops over all neighbors one more time to examine the effect of their displacement along the 3 coordinate axes one at a time. Similar to the CPU-based implementation, whenever the pair of force components ±δF cav = ±4πγ i (R off i ) 2 /(N δr) need to be added to the total solvation forces of two neighbor atoms a i and a j along the proper coordinate axis, the write operation that modifies F cav i of a i happens atomically in the shared memory. This ensures mutual exclusion between threads of the same block. On the other hand, the write operation that modifies F cav j of a j happens atomically in the global memory to ensure mutual exclusion between blocks of the same grid. Figure 7 shows the multi-threading scheme for the GPU-parallel algorithm. The algorithm is implemented using NVIDIA's compute-unified device architecture (CUDA). Kernel invocation is carried out synchronously within the default CUDA stream, hence synchronization between blocks is automatically guaranteed, while barrier synchronization is needed between threads of the same block.
GPU Optimization. The optimization attempts can be categorized as memory, execution, instruction, and flow-control optimization.
• Memory optimization is the most effective of all, as demonstrated by the results in the Section 5.2. In contrast to the CPU-parallel algorithm that makes most references through the cached global memory, the GPUparallel algorithm transfers the coordinates, radii, solvation parameters, and neighbor index lists for each atom into the shared memory to minimize the number of global memory references. The variables that are exclusive to the threads, on the other hand, such as the exposure states or sample coordinates are allocated in the registers. However, the limited amount of shared memory and register resources on the streaming multiprocessor (SM) imposes a restriction on the number of resident blocks on the SM and can adversely affect thread occupancy at any time during the simulation. Therefore, one needs to avoid excessive variable definitions within the scope of the GPU kernels.
• For execution level optimization, the kernels should be executed with proper granularity to maximize SM thread occupancy. Specifying a larger number of threads per block generally contributes to latency hiding, but is limited by the architecture as well as the on-chip memory resources. The number of threads is the same as the sample size N = |Q| a proper choice of which is a trade-off between accuracy and performance.
• For instruction level optimization, the transcendental math functions are converted to their intrinsic alternatives that are executed on the special function units (SFU) of the CUDA cores.
• Flow-control optimization is realized by avoiding multiple execution paths within the same block, which might lead to thread divergence and serialization within the same warp. In particular, when checking for overlaps between neighbor atoms and sample points, the conditional (e.g., if/else/then) statements are set in such a way that one of the two execution paths is always null.
The near-optimal conditions are reached by successive experimentation and modification of the code. For more information regarding the GPU architecture and terminology, see Appendix C.2.

Results & Discussion
This section presents a preliminary assessment of the model and implementation enhancements from Protofold I [46][47][48] to Protofold II. The folding process is simulated and assessed at multiple levels, ranging from the formation of secondary structural elements (e.g., α−helix coiling or β−strands formation) from an open chain to tertiary interactions between secondary elements or across larger domains that can be assumed to be rigid in real protein examples.  In Section 5.1 we discuss the impact of introducing solvation effects on the folding process of secondary structural elements starting from different initial conditions. We present some performance measures in Section 5.2 to validate the practical benefits of algorithmic improvements (e.g., coordinate hashing) as well as implementation improvements (i.e., CPU-and GPU-parallel computing). Finally, we look at a few real protein molecules in Section 5.3 and examine the energy variations in the neighborhood of the native structures.

The Folding Process
The simplest structural elements that are ubiquitous across many protein domains are α−helices and β−strands. Here we start by considering simple test runs on relatively short (e.g., 10-20 residues long) peptide chains made of Ala residues (typically used as the benchmark AA type in its most common L-stereoisomeric form) to visualize their compliance into secondary structural elements.
Alpha Helix Formation. First, we run four tests on a 15residue chain starting from two different initial conditions, for both of which we simulate the folding process without and with solvation effects taken into account: 1. Starting from the (slightly pre-coiled) initial conditions φ 0 i = ψ 0 i = +10 • for all 1 ≤ i ≤ 15 the chain folds into a left-handed α−helix as depicted in Fig. 8. The energy variation during KCM iterations is given in Fig. 10. 2. Starting from the (slightly pre-coiled) initial conditions φ 0 i = ψ 0 i = −10 • for all 1 ≤ i ≤ 15 the chain folds into a right-handed α−helix as depicted in Fig. 9. The energy variation during KCM iterations is given in Fig. 11.
The folding process in vacuum, i.e., without considering the solvation effects, emulates the behavior in the absence of the polar solvent, e.g., in membrane proteins extended along the nonpolar lipid bilayer or in secondary structural elements wrapped inside the hydrophobic core of globular proteins [3]. On the other hand, the folding process in water, i.e., with the presence of the solvation effects in addition to the intramolecular interactions, emulates the formation of elements that reside at the hydrophilic surface of globular proteins [3].
In the case of α−helix formations in Figs. 8 and 9, the hand of the initial coil determines the hand of the final helix. 11 This is due to the gradient descent nature of the KCM search algorithm that tends to converge to different local minima depending on the initial state. The effects of the solvation are hard to observe in these examples with the energetically unchallenged helical structures due to proper stacking of the atoms favored by all considered effects. In both cases the van der Waals and solvation effects work in the same direction until the steric clash prevents the helix to coil further. Figures 10 and 11 are plots of the free energy variations ver- sus KCM iteration number for the four runs described above.
Note that in all four cases (top and bottom plots) the solvation energy is evaluated and plotted, but only in two of them (bottom plots) its effects are applied to deform the chain.
For both left and right-handed helix formation, the inclusion of solvation effects clearly changes the folding pathway and increases the number of iterations before convergence from around 150 to 300. However, in either case the solvation free energy changes are not as significant as those of intramolecular (particularly van der Waals) effects. Another important observation is that the right-handed α−helix exhibits a notably more stable conformation than the left-handed α−helix with about ∼ 40−50 kcal per mol lower total free energy state-to be accurate, 43.8 and 45.9 kcal per mol for the entire chain, i.e., 2.9 and 3.1 kcal per mol per AA residue, without and with solvation effects, respectively. Although this is qualitatively consistent with the expectation of righthanded coiling being favored by L-alanine chains, the energy differences are higher than the ones reported in earlier studies (e.g., MD results in [86]). However, a meaningful comparison would require using identical simulation parameters, which is beyond the scope of this paper. The final dihedral angles for all 15 Ala residues corresponding to the folded (i.e., stable) conformations, obtained after a large enough number of iterations, are given in Table 4.
Ramachandran Plots. To examine the local effects of energetics, Ramachandran plots (for a pair of Ala residues in tandem) are generated by Protofold II using the energy-field presented in Section 2.2. The plots in Fig. 12 show the energy variations across different pairs (φ, ψ) of dihedral angles, without and with considering the solvation effects. The vertical high-energy regions in the middle corresponds to prohibitive steric clashes between the atoms. The low-energy regions around it, on the other hand, correspond to the ge- ometric relations between consecutive peptide planes that, when repeated for a segment of multiple residues along the chain, create secondary structural elements such as coiled α−helices and flat β−strands. Although the solvation effects do not significantly alter the shape of the energy profile, there is a certain amount of noise added due to the discrete nature of the enumeration algorithm presented in Section 3.4.
To observe the effects of solvation, one needs to carry out more extensive simulations on larger data sets with different chain lengths and various initial conditions. We carried out KCM runs on 2,000 independent polyalanine chains of random lengths in the range 10 ≤ m ≤ 20 starting from random initial angles in the range −90 The tests were run separately without and with considering the solvation effects using the same random seed. The resulting dihedral angles for all 27,717 Ala residues of all chains 12 are plotted in Fig. 13. One can observe multiple concentration areas that clearly correspond to left-and right-handed α−helices and flat β−strands, the former two helical folds being more populated on the plots. Zooming further on the two α−regions reveals multiple local minima where the points are concentrated more, which correspond to different subtypes of α−helices. Comparing the two plots in Fig. 13 reveals that the solvation effects do not significantly change the locations of the local minima. However, the energy profile is relaxed around the local minima where it has sharp cracks traced by the concentrated points along the valleys of the intramolecular energy landscape.
In water In vacuum In water

Computation Times
Next, we demonstrate the substantial performance improvements in Protofold II as a result of introducing the algorithms and data structures presented in Section 3. The running times are reported and compared on two computer systems, namely: • C-1: Dell Precision T7500 workstation with an Intel Xeon E5645 CPU (12 cores, clock rate 2.40 GHz, and host memory 24 GB). The system is equipped with one NVIDIA Quadro 4000 GPU (256 CUDA cores with compute capability (CC) 2.0 and device memory 2 GB).
• C-2: Dell Precision T7600 workstation with an Intel Xeon E5-2687W CPU (32 cores, clock rate 3.10 GHz, and host memory 64 GB). The system is equipped with two graphics cards: a NVIDIA Quadro K5000 GPU (1,536 CUDA cores with CC 3.0 and device memory 4 GB) and a NVIDIA Tesla K20C GPU (2,496 CUDA cores with CC 3.5 and device memory 5 GB).
Effect of Hashing. Figure 14 shows the running times of the force computation step (on C-1) in a single KCM iteration for folding polyalanine chains of different lengths with and without 3D hashing presented in Seciton 3.2 both in vacuum (i.e., without considering solvation effects) and in water (i.e., with considering solvation effects). In both cases, the results show a significant reduction in the running times with hashing (e.g., up to 4.6× in vacuum and 1.5× in water for m = 60 Ala residues), and the difference scales with the size of the molecule. Nevertheless, the solvation force computations remain the bottleneck of the simulation (using sequential CPU implementation) and adversely affect the speed-up gained from hashing.
Parallel Computing. The first 3 columns on the left in Fig. 15 show the sequential CPU running times (on C-1) for the electrostatic, van der Waals, and nonpolar solvation Figure 14: The effect of 3D hashing on the force computation times of a 60-residue polyalanine chain with and without solvation effects (on C-1). Figure 15: Sequential running times for electrostatic, van der Waals, and solvation forces, and CPU-parallel running times of the latter for a 1,200-residue polypeptide chain (on C-1).
forces in a single KCM iteration. The results were obtained for a polypeptide chain composed of 1, 200 residues that are randomly selected from Ala, Cys, and Ser AAs-the latter two being equivalent to replacing the H in the nonpolar Ala side chain with SH and OH, respectively, resulting in polar side chains. It is clearly observed that with a single CPU core, the first two terms take a very small fraction of the total time (around 5−6% of the total force computation time per iteration) and the solvent effects are clearly the bottleneck. The same column chart in Fig. 15 also shows the running times and corresponding speed-ups (on C-1) for the CPU-parallel computation of the solvation force using up to 12 CPU cores. An almost linear speed-up is achieved by increasing the number of processors (e.g., ∼ 10× with 12 cores). However, the solvation force calculation is still about an order of magnitude slower than that of the other two force types. Figure 16 compares the running times and corresponding speed-ups (on C-1) of the CPU-and GPU-parallel implementations for polypeptides of different lengths, ranging from m = 200 (i.e., ∼ 2K atoms) to 1, 200 AAs (i.e., ∼ 13K atoms). To depict the importance of memory optimization presented in Section 4.3, the running times are shown for two different cases; namely, one that uses global memory for communications between all threads of all blocks, and the optimized code making extensive use of shared memories for communications between threads within the same block. It is interesting to note that when the GPU shared memory is not utilized, the results do not show an improvement over the 12-core CPU implementation due to large global memory access latencies. However, proper shared memory usage results in a huge performance improvement that scales by protein size, e.g., from 20× for m = 200 AAs to 70× for m = 1, 200 AAs with respect to the sequential run on a single CPU.
The computations are repeated for even larger molecules in Fig. 17, ranging from m = 1, 000 AAs (i.e., ∼ 10K atoms) to 6, 000 AAs (i.e., ∼ 64K atoms), enabled by leveraging a more powerful machine (i.e., C-2). For the case of m = 1, 000 AAs, C-2 yields a two-fold speed-up on the CPU and a three-fold speed-up on the GPU compared to C-1. As depicted in Fig. 17, significantly higher and more consistent CPU speed-ups of 16-18× and GPU speed-ups of 90-100× (for Quadro K5000) and 270-290× (for Tesla K20C) are observed. These observations imply proper scalability of the data-parallel implementation with molecular size.
Even for molecules with tens of thousands of atoms, each force-field computation takes only a fraction of a second per iteration. This enables fast KCM simulation of folding for large proteins over extended periods of time via Protofold II, which wouldn't be tractable via Protofold I [46][47][48].

Real Examples
Having considered the folding of secondary structural elements with complete flexibility (i.e., each peptide plane being treated as a separate rigid body), we proceed to study the tertiary interaction between larger rigid units in real proteins. We report on the following case studies:  and Troponin-C (PDB: 2JNF).
• The interactions between multiple rigid domains that are connected via flexible loops within the containing monomeric unit are examined. The example of Gamma-B Crystallin (PDB: 1GCS) is used for this purpose.
These PDB structures that are obtained from X-ray crystallography (e.g., 1TES and 1GCS) do not contain the H atoms, hence are first preprocessed using Duke University's MolProbity server [87,88] which predicts and adds the H atoms to the coordinates information before importing the structure into Protofold II. The PDB structures that are obtained from nuclear magnetic resonance (NMR) spectrometry (e.g., 2JNF), on the other hand, already contain the H atoms positions. The size of the molecules, i.e., the number of AA residues m and the number of atoms n and n with and without the H atoms included, respectively, are given in Table 5.
Secondary Structural Interactions. Let us first consider the energy variations when an α−helix of an α−domain    is reoriented from its stable conformation with respect to the rest of the bundle, as illustrated in Figs. 18 and 20. Myoglobin (PDB: 1TES) is an oxygen binding muscle protein that is composed of a single 'globin fold' domain, which is an α−domain motif consisting of a bag of 8 α−helices per domain (denoted A through H) arranged at ∼ +90 • and +50 • angles with respect to each other, as shown in Fig.  18. This arrangement creates a hydrophobic pocket in the interior that wraps the stabilizer co-factor known as 'heme group' [89]. Assuming that the helices are rigid, we examine the energy variations due to dihedral rotations at the loops that connect the two end α−helices; namely, local changes in (φ i , ψ i ) for i = 21 and i = 125 where the A and H helices are hinged, respectively.
Troponin-C (PDB: 2JNF) is a calcium binding muscle protein that is composed of two 'EF-hand' domains, which are α−domain motifs consisting of a bundle of 4 α−helices per domain (denoted A through H) arranged in an up-and-down anti-parallel conformation [89], as shown in Fig. 20. There is a long pseudo-helical segment that connects the two globular domains. Assuming that the helices are rigid, we examine the energy variations due to dihedral rotations at the loops that connect the two end α−helices; namely, local changes in (φ i , ψ i ) for i = 24 and i = 145 where the A and H helices are hinged, respectively. Figures 19 and 21 show the free energy variations in vacuum (i.e., without considering solvation effects) and in water (i.e., with considering solvation effects) for the two protein domains as the end α−helices A and H of each are rotated within a range of ±10 • of the native (φ i , ψ i ) at the hinged loops. 13 In all four cases the energy model of Protofold II exhibits a local minima within ±1 • of the native conformation. We also observe that the solvation effects contribute a significant amount to the total free energy; however, the SASA variations are so small in the considered neighborhood that the location of the local minima is almost unchanged. The van der Waals effects appear to be dominant in this neighborhood, manifested as the shape complementarity between the ridges and grooves of the mobile α−helix and those of the other helices in the bundle [90]. However, when variations across larger angular ranges are considered, the solvation effects are expected to play a more determining role.
Tertiary Structural Interactions. Lastly, we consider the energy variations when a rigid domains of a protein is reoriented with respect to another domain against which it is packed into a stable structure.
Gamma-B Crystallin (PDB: 1GCS) is an eye-lens protein that is made of two similar domains that are 40% identical in sequence [89]. Each domain is composed of two anti-parallel β−sheets each made of 4 β−strands with the same arrangement topology [89], as shown in Fig. 22. Assuming that these domains are rigid, we analyze the energy variations due to rotating one domain with respect to the other at one of the 13 The colormaps are generated using a nonlinear but consistent contouring scheme. residues that belong to the connecting loop (e.g., i = 81). To observe the global effects of solvation, we allow both dihedral angles (φ 81 , ψ 81 ) to vary over the entire range of ±180 • from the native values. To facilitate visualization, this time we consider only one angle's variation at a time. Figures 23 and 24 show the variations of the different energy terms. 14 It appears that the van der Waals effects are dominant in determining the profile of the energy well near the native conformation, which can be attributed to the extensive contact interface between the two domains in Fig. 22. The electrostatic and solvation effects substantially change the energy landscape thus can dramatically affect the folding pathway. However, they do not cause a significant change in the location of the energy minimum. Once again, the solvation energy is noisier due to the discrete nature of the enumeration algorithm in Section 3.4. Another interesting observation is that the electrostatic energy has discontinuities due to the cut-off approximation when pairs of atoms are farther than 9.0Å. Although it does not seem to affect the minimum location in this example, larger cut-off distances might be necessary for analyzing large proteins, since the accumulation of the pairwise errors grows quadratically with the number of atoms.

Conclusion
The KCM approach to protein folding [44,45] originally implemented into Protofold I [46][47][48] provides a promising fast alternative to the popular MD simulation and MC sampling methods by 1. modeling the protein chain as a kinematic linkage with restricted DOF to which the well-studied principles from mechanism synthesis and robotics can be readily applied; and 2. replacing the 2nd-order dynamic response with 1st-order kinetostatic integration of the equations of motion to facilitate convergence to the free energy minima.
In the present work, we introduced major model and implementation improvements in Protofold II by 14 To enabel logarithmic plots, each energy ordinate is offset by a constant value to shift its minimum to (the arbitrary positive value of) 100 kcal per mol.  1. incorporating the solvation effects that characterize the hydrophobic effect, i.e., the entropic changes due to cavity formation in the aqueous solvent; 2. taking advantage from efficient auxiliary algorithms and data structures to improve the computational complexity from O(n 2 ) to expected O(n); and 3. implementing fast and relatively accurate evaluation of the SASA and its gradient for solvation energy-and force-field computation, respectively, in parallel on both CPU and GPU.
The presented enumeration algorithm for the latter provides a fast approximate method, in which the degree of accuracy is traded off with the performance by a proper choice of the sample size. We argued that the inclusion of the solvation free energy into the mix can significantly affect the folding pathway for water-soluble proteins whose folding in vivo is dominated by such effects. We also demonstrated that the performance gain of the GPU-accelerated implementation scales properly with the number of atoms, achieving up to two orders of magnitude in speed-ups after memory optimization.
Protofold II has been completely rearchitectured and is evolving into a versatile analysis toolbox for studying the kinematic and structural behavior of molecular chains in protein engineering applications such as the design of nanomanipulators [53,54] among the ongoing projects. A hybrid force-field model was used, composed of 1. the AMBER model [67] (for noncovalent interactions); and 2. supplemental terms similar to the CHARMM model [78] and the GROMOS model [64] (for solvation effects) except for the probabilistic SASA estimation [63] replaced with our own surface sampling algorithm.
This model is by no means versatile enough to enable addressing the ultimate goal of arriving from sequence to 3D structure at a click of a button. Even though predicting folding of real 3D structures requires further developments, this study represents a major step toward this goal.

B Prefix Computation
The prefix sum problem is fundamental to numerous important algorithms, and is defined as follows. Given a finite ordered sequence of elements X = (x 1 , x 2 , · · · , x n ) ∈ Σ n and an arbitrary binary operator ⊕ : Σ × Σ → Σ that is O(1)−time computable and associative (i.e., (x ⊕ y) ⊕ z = x ⊕ (y ⊕ z)), compute another sequence Y = (y 1 , y 2 , · · · , y n ) ∈ Σ n where y 1 = x 1 , y 2 = x 1 ⊕ x 2 , · · · , y n = x 1 ⊕ x 2 ⊕ · · · ⊕ x n ; in other words y i = y i−1 ⊕ x n for 1 ≤ i ≤ n where y 0 is the left-identity element (i.e., y 0 ⊕ x = x for all x ∈ Σ). It is trivial to show that the prefix sums can be computed sequentially in O(n), which is optimal. In addition, there are work-optimal parallel algorithms with a total computational work of T P = O(n) that carry out the prefix computation in T = O(log n) time using P = O(n/ log n) CREW PRAM or in T = O(log n/ log log n) time using P = O(n log log n/ log n) CRCW PRAM processors with common conflict resolution [70]-see Appendix C.1 for details regarding PRAM.

C.1 Abstract Machines
The parallel random-access machine (PRAM) is a sharedmemory abstract parallel computation model, typically assigned with the exclusive/concurrent-read (ER/CR) and exclusive/concurrent-write (EW/CW) attributes [85]. The most common attributes are CREW and CRCW, noting that multiple processors can concurrently read a memory cell but only one can write at a time to prevent race conditions. To enable concurrent writing, one needs to resolve possible conflicts with typical mechanisms such as 1) 'common' meaning that all processors attempt to write the same value); 2) 'arbitrary' meaning that one processor's attempt succeeds at random; and 3) 'priority' meaning that the processors are prioritized by a prespecified order. Note that a CREW algorithm can always run in the same (if not fewer) number of steps on a CRCW machine.

C.2 GPU SIMT Model
A typical CUDA GPU program proceeds as follows. The data is first transferred from the CPU (i.e., host) memory to the GPU (i.e., device) memory. The host application invokes the so-called kernels on the GPU with specified granularity, i.e., issuing a 1D, 2D, or 3D grid of blocks, each block being a 1D, 2D, or 3D array of threads that are sent in groups of 16 or 32 (called 'warps') to one of the streaming multiprocessors (SM). The threads within the same block can access the fast shared memory banks on the SM, and communication across blocks is done using the global memory. The computed results are transferred from the device memory back to the host memory.
There are different types of GPU memory locations, classified into two groups: 1) device (i.e., off-chip) memory in-cluding global and local memories; and 2) on-chip memory including shared memory, cache, and registers. The access latencies to the on-chip are much less (around 100× faster) than those of the off-chip memory.