the objective of the present study is 
 to simulate the behavior of mesoscopic systems based on an all - atom 
 formulation at which the basic physics is presumed known . traditional 
 molecular dynamics ( md ) is ideal for such an approach if the number 
 of atoms and the time scales of interest are limited . however , ribosomes , viruses , mitochondria , and nanocapsules for 
 the delivery of therapeutic agents are but a few examples of mesoscopic 
 systems that can provide a challenge for conventional md . in this 
 paper , we develop a physics - based algorithm that accounts for interactions 
 at the atomic scale and yet makes accurate and rapid simulations for 
 supramillion atom systems over long time scales possible . typical 
 coarse - graining ( cg ) methods include deductive multiscale analysis 
 ( dma ) , inverse monte carlo , boltzmann inversion , elastic network 
 models , or other bead - based models . dma methods derived from the n - atom liouville equation 
 ( le ) show great promise in achieving accurate and efficient all - atom 
 simulation . the main theme of that work was to construct and exploit the multiscale 
 structure of the n - atom probability density (,t ) for the positions and momenta of the n atoms ( denoted  , collectively ) as it evolves over time t. most of the analysis focused on friction dominated , noninertial 
 regime , which is considered here as well . however , in these methods 
 ensembles of all - atom configurations were required for evolving the 
 cg variables . the approach introduced here avoids the need to construct 
 these ensembles by coevolving the all - atom and cg states in a consistent 
 way , and in the spirit of dma - based methods , it does not make any 
 conjectures on the form of the cg dynamical equations and the associated 
 uncertainty in the form of the equations . a main theme of the present 
 approach is the importance of coevolving the cg and microscopic states . this feature distinguishes our method from others which , for example , 
 require the construction of a potential mean force using ensembles of micostates ; a challenge for such methods is that 
 the relevant ensembles are not known a priori since they are controlled 
 by the cg state whose evolution is unknown and is in fact the objective 
 of a dynamics simulation . other multiscale methods , built on the projection 
 operator formalism , require the construction of memory kernels . this 
 is typically achieved via a perturbation approach to overcome the 
 complexity of the appearance of the projection operator in the memory 
 kernels . a first step in the present approach is the introduction of a set 
 of cg variables  related to  via  = ( ) 
 for specified function ( ) . when this dependence 
 is well chosen , the cg variables evolve much more slowly than the 
 fluctuations of small subsets of atoms . with these cg variables , the n - atom le was solved perturbatively in terms of  , the ratio of the characteristic time of the fluctuations of small 
 clusters of atoms , to the characteristic time of cg variable evolution . this is achieved starting with the ansatz that  depends on 
  both directly and , via  , indirectly . the theory 
 proceeds by constructing ( , ;t ) perturbatively in  , i.e. , by working in the space of functions 
 of 6n + ncg variables ( where ncg is the number of variables 
 in the set  ) . to advance the multiscale approach , we here introduce 
 trotter factorization into the analysis . through trotter factorization , 
 the long - time evolution of the system separates into alternating phases 
 of all - atom simulations and cg variable updating . efficiency of the 
 method follows from a hypothesis that the momenta conjugate to the 
 cg variables can be represented as a stationary process whose statistics 
 are preserved in an interval of time short to the characteristic time 
 of cg dynamics . the net result is a computational algorithm with some 
 of the character of our earlier md / opx method but with better control on accuracy , higher efficiency , and more 
 rigorous theoretical basis . here we develop the algorithm and discuss 
 its implementation as a computational platform , discuss selected results , 
 and make concluding remarks . the newtonian description of an n - atom system is provided by the 6n atomic positions 
 and momenta , denoted  collectively . the phenomena of interest 
 involve overall transformations of an n - atom system . 
 while  contains all the information needed to solve the problem 
 in principle , here it is found convenient to also introduce a set 
 of cg variables  , that are used to track the large spatial 
 scale , long time degrees of freedom . for example ,  could describe 
 the overall position , size , shape , and orientation of a nanoparticle . 
 by construction , a change in  involves the coherent deformation 
 of the n - atom system , which implies that the rate 
 of change in  is expected to be slow . this slowness implies the separation of time scales that provides 
 a highly efficient and accurate algorithm for simulating n - atom systems . with this unfolded description (  ,  ) , the newtonian dynamics takes the form12for unfolded 
 liouvillian , such that34here k is the cg velocity associated with the kth 
 cg variable . equations 1 and 2 have the formal solution5for initial data indicated 
 by subscript o and evolution operator s(t ) = . by taking the unfolded liouvillian , the time 
 operator now takes the form6since   and   do 
 not commute , s(t ) can not be factorized 
 into a product of exponential functions . however , trotter s 
 theorem ( also known as the lie product 
 formula ) can be used to factorize the 
 evolution operator as follows:7by setting t / m to be equal to the discrete time step  , the stepwise 
 operator becomes8let the stepwise operators smicro and smeso correspond to   and , respectively . then s(n ) takes the form9by replacing smicro(/2 ) by smicro()smicro(/2 ) to 
 the right - hand side , eq 9 becomes10since we are interested in the long - time evolution 
 of a mesoscopic system , we can neglect the far left and right end 
 terms , smicro(/2 ) and smicro(/2 ) , respectively , to a good approximation . therefore , we can define the stepwise time operator as11 in the next section , we show how this factorization implies a computational 
 algorithm for solving the dynamical equations for  and . a key to the efficiency of the mutiscale trotter factorization 
 ( mtf ) method is the postulate that the momenta conjugate to the cg 
 variables can be represented by a stationary process over a period 
 of time much shorter than the time scale characteristic of cg evolution . thus , in a time period significantly shorter than the increment  
 of the stepwise evolution , the system visits a representative ensemble 
 of configurations consistent with the slowly evolving cg state . this 
 enables one to use an md simulation for the microscopic phase of the 
 stepwise evolution that is much shorter than  to integrate 
 the cg state to the next cg time step . for each of a set of time intervals much less than  , the friction dominated system experiences 
 the same ensemble of conjugate momentum fluctuations . thus , if  
 is the time for which the conjugate momentum undergoes a representative 
 sample of values ( i.e. , is described by the stationarity hypothesis ) , 
 then the computational advantage over conventional md is expected 
 to be /. the two - phase updating for each time - step 
  was achieved as follows . for the smicro( ) phase , this yields a time - series 
 for  and hence . for all systems simulated here ,  
 was found to be a stationary process ( see figure 10 ) . therefore , md need only be carried out for a fraction of 
  , denoted . this and the slowness of the cg variables 
 are the source of computational efficiency of our algorithm . for the smeso phase updating in the friction dominated 
 regime , the  time series constructed in the microphase is used 
 to advance  in time as follows12due to stationarity , the integral on the right - hand side reduces 
 to /tt+dt (t ) ( see figure 10 ) . the expression for  depends on the choice 
 of cg variables . in this work , we used the space - warping method that maps a set of atomic coordinates to a set of cg variables that 
 capture the coherent deformation of a molecular system in space . in 
 the space - warping method , the mathematical relation between the cg 
 variables and the atomic coordinates is13here k is a triplet of indices , i is the atomic index , ri is the cartesian position vector for atom i , and k is a cartesian 
 vector for cg variable k. the basis functions uk are constructed in two 
 stages . in the first stage , they are computed from a product of three 
 legendre polynomials of order k1 , k2 , and k3 for the x , y , and z dependence . 
 in the second stage , for instance , the zeroth order 
 polynomial is u000 , the first order polynomial 
 forms a set of three basis functions : u001 , u010 , u100 , and 
 so on . furthermore , the basis functions depend on a reference configuration r which is updated periodically ( once every 10 
 cg time steps ) to control accuracy . the vector si represents the atomic - scale 
 corrections to the coherent deformations generated by k. introducing cg variables this 
 way facilitates the construction of microstates consistent with the 
 cg state . this is achieved by minimizing mii2 with respect to k. the result is that the cg variables 
 are generalized centers of mass , specifically14with mi being the mass 
 of atom i. for the lowest order cg variable , u000 = 1 , which implies 000 is the center of mass . as the order of the polynomial increases , the cg variables capture more information from the atomic scale , but 
 they vary less slowly with time . therefore , the space warping cg variables 
 are classified into low - order and high - order variables . the former 
 characterize the larger scale disturbances , while the latter capture 
 short - scale ones . equation 14 implies that k = ukipi / miuki2 , where pi is a vector of momenta for 
 the i atom . with (t +  ) computed via eq 12 , the two - phase  update is completed and 
 this cycle is repeated for a finite number of discrete time steps . details on the necessary energy minimization and equilibriation needed 
 for every cg step was covered in earlier work . this two - phase coevolution algorithm was 
 implemented using namd for the smicro phase within the framework of the dms 
 software package . numerical computations were performed with the aid of loos , a lightweight object - oriented structure library 
 for md simulation analysis . all simulations 
 were done in vacuum under nvt conditions to assess 
 the efficiency and accuracy of the algorithm . this iron binding 
 protein is composed of a distal and two proximal lobes ( shown in figure 1 ) . two free energy minimizing conformations have 
 been demonstrated experimentally : diferric with closed proximal lobes 
 ( pdb code 1lfg ) and apo with open ones ( pdb code 1lfh ) . here , we start 
 with an open lactoferrin structure and simulate its closing in vacuum . a total of 10  3 cg variables are used ( i.e. , k1 + k2 + k3  2 ) to coarse - grain the protein , and the highest - order 
 basis functions are therefore quadratic . this level of coarse - graining 
 captures translation , rotation , contraction , and bending of a macromolecule . the rmsd for lactoferrin is plotted as a function 
 of time in figure 4 ; it shows that the protein 
 reaches equilibrium in about 5 ns . this transition leads to a decrease 
 in the radius of gyration of the protein by approximately 0.2 nm as 
 shown in figure 5 . snapshots of lactoferrin 
 protein in its open state at t = 0 ns ( a ) and its 
 closed state at t = 19.6 ns ( b ) . the second system simulated is a triangular structure of 
 the nudaurelia capensis omega virus ( nv ) capsid protein ( pdb code 1ohf ) containing three protomers ( shown in 
 figure 2 ) . starting from a deprotonated state 
 ( at low ph ) , the system was equilibriated using an implicit solvent . a total of 10  3 cg variables are used ( i.e. , k1 + k2 + k3  2 ) to coarse - grain the structure . as 
 a result , the proteins shrink in vacuum after a short period of equilibriation . snapshots of nv triangular structure before ( initial 
 state at t = 0 ns ) ( a ) and after ( shrinking at t = 3.0 ns ) ( b ) contraction due to strong protein the third system simulated 
 is an assembly of l1 proteins for the human papillomavirus ( hpv ) . a multiscale simulation is initiated by configuring 
 a set of nine pentamers ( in contrast to the twelve pentamers constituting 
 the t = 1 structure ) . partial structures are of interest 
 as they arise along the self - assembly pathway and their stability 
 could provide a kinetic road block for their complete self - assembly . the space - warping variables were used to coarse - grain the system at 
 the pentamer level , i.e. , each of the nine pentamers was replaced 
 by its center of mass . this level of coarse - graining captures the 
 translational movement of the pentamers as they move closer to one 
 another in order to minmize their free energy . the number of hydrogen 
 bonds increases in time , which is captured by the mtf algorithm ( shown 
 in figure 7 ) . for the first and second pentamers , 
 the z - component of their centers of mass and the 
 autocorrelation functions of their y - components are 
 shown in figures 8 and 9 , respectively . human papillomavirus ( hpv ) capsid - like structure of l1 
 proteins ( left ) is divided into nine pentameric subsystems ( right ) , 
 each of which has its own set of cg variables ( taken here to be the 
 center of mass of each pentamer , i.e. , 000 ) . rmsd variation as a function of time for a series 
 of three md and one mtf runs . temporal evolution of the radius of gyration 
 of nv computed using md and mtf . time variation in the number of hydrogen bonds obtained via md and 
 mtf for the hpv t = 1 assembly of l1 proteins . z - component of the centers 
 of mass for the first ( left ) and second ( right ) pentamers increases 
 in time as the nine pentamers in the hpv assembly move closer to one 
 another . time decay in the autocorrelation function ( acf ) 
 for the y component of the center mass for the first 
 ( left ) and second ( right ) pentamer in the hpv assembly of proteins . on the basis of the convergence 
 of the time integral of  ( shown in figure 10 ) , the smicro phase was chosen to consist of 5 ps for lfg , 10 
 ps for nv , and 1 ps for hpv . the cg time step ,  , was 
 taken to be 12.5 ps for lfg , 25 ps for nv , and 10 ps for hpv . this yielded a speedup ( with respect to md ) of 1.32 , 2.15 , and 7.15 , 
 respectively . evidence for the validity of the stationarity hypothesis 
 shown via the convergence of ( 1/)0(t ) dt as a function of  for cg variables selected 
 from among those used in simulating the contraction of nv . initially the integral experiences large fluctuations because with 
 small  , only a relatively few configurations are included in 
 the time average constituting the integral , but as  increases , 
 the statistics improves , and the integral becomes increasingly flat . ( b ) ( a ) time integral of  for a low - order cg 001 . mesoscopic 
 systems express behaviors stemming from atom  atom interactions 
 across many scales in space and time . earlier approaches based on 
 langevin equations for coarse - grained variables showed agreement with 
 md without comprimising accuracy and captured key atomic scale details . however , such an approach requires the construction of diffusion 
 factors , a task that consumes significant computational resources . this is because of the need to use large ensembles and construct correlation 
 functions . the multiscale factorization method used here introduces 
 the benefits of multiscale theory of the le . here we revisit the multiscale 
 method within our the context of mtf which avoids the need for the 
 resource - consuming construction of diffusion factors , and thermal 
 average and random forces . the cg variables for the mesoscopic systems 
 of interest do have a degree of stochastic behavior . in the present 
 formulation , this stochasticity is accounted for via a series of md 
 steps used in the phase of the multiscale factorization algorithm 
 wherein the state of the system is evolved via , i.e. at constant value of the cg variables . the mtf algorithm can be further optimized to produce greater speedup 
 factors . in particular , the results obtained here can be significantly 
 improved with the following : ( 1 ) after updating the cgs in the two - phase 
 coevolution trotter cycle , one must fine - grain , i.e. develop , the 
 atomistic configuration to be used as an input to md . we have shown 
 that the cpu time to achieve this fine - graining can be dramatically 
 reduced via a constraint method that eliminates bond length and angle 
 strains , ( 2 ) information from earlier steps in discrete time coevolution 
 can be used to increase the time step and achieve greater numerical 
 stability ; while this was demonstrated for one multiscale algorithm , it can also be applied to mtf , and ( 3 ) the time 
 stepping algorithm used in this work is the analogue of the euler 
 method for differential equations , and greater numerical stability 
 and efficiency could be achieved for a system of stiff differential 
 equations using implicit and semi - implicit schemes .
mesoscopic n - atom systems derive their structural and dynamical properties 
 from processes coupled across multiple scales in space and time . a 
 multiscale method for simulating these systems in the friction dominated 
 regime from the underlying n - atom formulation 
 is 
 presented . 
 the method integrates notions of multiscale analysis , trotter 
 factorization , and a hypothesis that the momenta conjugate to coarse - grained 
 variables constitute a stationary process on the time scale of coarse - grained 
 dynamics . 
 the method is demonstrated for lactoferrin , nudaurelia capensis 
 omega virus , and human papillomavirus to assess its accuracy .