Simulation Of Large Systems Research Articles

Advances in Atomic-Level Simulations of Large-Scale Functional Motions of Membrane Transporters

Membrane transporters are specialized molecular devices that use various forms of cellular energy to drive active transport of their specific substrates across the membrane. Their fundamental role in diverse key biological processes has placed them among central drug targets, furthering widespread interest in their biophysical and mechanistic studies at a molecular level. Large-scale conformational changes are central to the function of membrane transporters. Description of these structural changes, however, requires sampling high-dimensional free energy landscapes that are inaccessible to conventional sampling techniques such as regular molecular dynamics (MD) simulations. We have recently developed a novel computational approach that, while being numerically expensive, has been the most efficient way to describe large-scale structural transitions in membrane transporters (as well as for any other macromolecular systems) using non-equilibrium methods employing system-specific collective variables, and a novel combination of several state-of-the-art sampling techniques [Moradi and Tajkhorshid, PNAS 110:18916-21 (2013); JCTC 10: 2866-80, 2014]. The approach is based on loosely coupled, multiple-copy MD simulations of large macromolecular systems preserving realistic representations of the systems in explicit membranes, and therefore relies on massive computing resources. Here we describe the application of the methodology to the study of several classes of membrane transporters, in order to characterize the inter-conversion of these molecular devices between the major conformational states necessary for their function, to characterize the free energy profiles associated with these transitions, and more importantly how chemical details such as ion/substrate binding drastically modulate the energy landscapes. The results of these simulations elucidate highly relevant mechanistic details of the function of membrane transporters providing a detailed structural basis for the experimentally observed phenomena.

Real-time simulation of large open quantum spin systems driven by dissipation

We consider a large quantum system with spins $\frac{1}{2}$ whose dynamics is driven entirely by measurements of the total spin of spin pairs. This gives rise to a dissipative coupling to the environment. When one averages over the measurement results, the corresponding real-time path integral does not suffer from a sign problem. Using an efficient cluster algorithm, we study the real-time evolution from an initial antiferromagnetic state of the two-dimensional Heisenberg model, which is driven to a disordered phase, not by a Hamiltonian, but by sporadic measurements or by continuous Lindblad evolution.

HONPAS: A linear scaling open‐source solution for large system simulations

HONPAS is an ab initio electronic structure program for linear scaling or O(N) first‐principles calculations of large and complex systems using standard norm‐conserving pseudopotentials, numerical atomic orbitals (NAOs) basis sets, and periodic boundary conditions. HONPAS is developed in the framework of the SIESTA methodology and focuses on the development and implementation of efficient O(N) algorithms for ab initio electronic structure calculations. The Heyd–Scuseria–Ernzerhof (HSE) screened hybrid density functional has been implemented using a NAO2GTO scheme to evaluate the electron repulsion integrals (ERIs) with NAOs. ERI screening techniques allow the HSE functional calculations to be very efficient and scale linearly. The density matrix purification algorithms have been implemented, and the PSUTC2 and SUTC2 methods have been developed to deal with spin unrestricted systems with or without predetermined spin multiplicity, respectively. After the self‐consistent field (SCF) process, additional O(N) post‐SCF calculations for frontier molecular orbitals and maximally localized Wannier functions are also developed and implemented. Finally, an O(N) method based on the density matrix perturbation theory has been proposed and implemented to treat electric field in solids. This article provides an overall introduction to capabilities of HONPAS and implementation details of different O(N) algorithms. © 2014 Wiley Periodicals, Inc.

Model reduction of descriptor systems using frequency limited Gramians

Model reduction is a process of approximating higher order original models by comparatively lower order models with reasonable accuracy in order to provide ease in design, modeling and simulation for large complex systems. Generally, model reduction techniques approximate the higher order systems for whole frequency range. However, some applications require approximation over a certain band of frequency than whole frequency range. Different model reduction techniques are available for descriptor systems. However, for limited frequency interval, no such work exists in the literature. A frequency limited balanced truncation method for general descriptor systems is proposed. The method is an extension of truncated balanced realization method for the general descriptor system. The proposed technique generalizes the results of Gawronski and Juang technique for large-scale descriptor systems using frequency interval Gramians. Simple algorithms are also given for preserving the stability of reduced-order models. The work also extends Poor Man׳s truncated balanced realization technique to include frequency limited Gramians for descriptor systems. Practical numerical examples are incorporated to show the successful application of the proposed method in the desired frequency range.

Efficient Maintenance and Update of Nonbonded Lists in Macromolecular Simulations.

Molecular mechanics and dynamics simulations use distance based cutoff approximations for faster computation of pairwise van der Waals and electrostatic energy terms. These approximations traditionally use a precalculated and periodically updated list of interacting atom pairs, known as the “nonbonded neighborhood lists” or nblists, in order to reduce the overhead of finding atom pairs that are within distance cutoff. The size of nblists grows linearly with the number of atoms in the system and superlinearly with the distance cutoff, and as a result, they require significant amount of memory for large molecular systems. The high space usage leads to poor cache performance, which slows computation for large distance cutoffs. Also, the high cost of updates means that one cannot afford to keep the data structure always synchronized with the configuration of the molecules when efficiency is at stake. We propose a dynamic octree data structure for implicit maintenance of nblists using space linear in the number of atoms but independent of the distance cutoff. The list can be updated very efficiently as the coordinates of atoms change during the simulation. Unlike explicit nblists, a single octree works for all distance cutoffs. In addition, octree is a cache-friendly data structure, and hence, it is less prone to cache miss slowdowns on modern memory hierarchies than nblists. Octrees use almost 2 orders of magnitude less memory, which is crucial for simulation of large systems, and while they are comparable in performance to nblists when the distance cutoff is small, they outperform nblists for larger systems and large cutoffs. Our tests show that octree implementation is approximately 1.5 times faster in practical use case scenarios as compared to nblists.

A Novel, Computationally Efficient Multipolar Model Employing Distributed Charges for Molecular Dynamics Simulations.

A truncated multipole expansion can be re-expressed exactly using an appropriate arrangement of point charges. This means that groups of point charges that are shifted away from nuclear coordinates can be used to achieve accurate electrostatics for molecular systems. We introduce a multipolar electrostatic model formulated in this way for use in computationally efficient multipolar molecular dynamics simulations with well-defined forces and energy conservation in NVE (constant number-volume-energy) simulations. A framework is introduced to distribute torques arising from multipole moments throughout a molecule, and a refined fitting approach is suggested to obtain atomic multipole moments that are optimized for accuracy and numerical stability in a force field context. The formulation of the charge model is outlined as it has been implemented into CHARMM, with application to test systems involving H2O and chlorobenzene. As well as ease of implementation and computational efficiency, the approach can be used to provide snapshots for multipolar QM/MM calculations in QM/MM-MD studies and easily combined with a standard point-charge force field to allow mixed multipolar/point charge simulations of large systems.

Ab initio nonadiabatic dynamics of multichromophore complexes: a scalable graphical-processing-unit-accelerated exciton framework.

Conspectus Although advances in computer hardware and algorithms tuned for novel computer architectures are leading to significant increases in the size and time scale for molecular simulations, it remains true that new methods and algorithms will be needed to address some of the problems in complex chemical systems, such as electrochemistry, excitation energy transport, proton transport, and condensed phase reactivity. Ideally, these new methods would exploit the strengths of emerging architectures. Fragment based approaches for electronic structure theory decompose the problem of solving the electronic Schrodinger equation into a series of much smaller problems. Because each of these smaller problems is largely independent, this strategy is particularly well-suited to parallel architectures. It appears that the most significant advances in computer architectures will be toward increased parallelism, and therefore fragment-based approaches are an ideal match to these trends. When the computational effort involved scales with the third (or higher) power of the molecular size, there is a large benefit to fragment-based approaches even on serial architectures. This is the case for many of the well-known methods for solving the electronic structure theory problem, especially when wave function-based approaches including electron correlation are considered. A major issue in fragment-based approaches is determining or improving their accuracy. Since the Achilles' heel of any such method lies in the approximations used to stitch the smaller problems back together (i.e., in the treatment of the cross-fragment interactions), it can often be important to ensure that the size of the smaller problems is "large enough." Thus, there are two frontiers that need to be extended in order to enable molecular simulations for large systems and long times: the strongly coupled problem of medium sized molecules (100-500 atoms) and the more weakly coupled problem of decomposing ("fragmenting") a molecular system and then stitching it back together. In this Account, we address both of these problems, the first by using graphical processing units (GPUs) and electronic structure algorithms tuned for these architectures and the second by using an exciton model as a framework in which to stitch together the solutions of the smaller problems. The multitiered parallel framework outlined here is aimed at nonadiabatic dynamics simulations on large supramolecular multichromophoric complexes in full atomistic detail. In this framework, the lowest tier of parallelism involves GPU-accelerated electronic structure theory calculations, for which we summarize recent progress in parallelizing the computation and use of electron repulsion integrals (ERIs), which are the major computational bottleneck in both density functional theory (DFT) and time-dependent density functional theory (TDDFT). The topmost tier of parallelism relies on a distributed memory framework, in which we build an exciton model that couples chromophoric units. Combining these multiple levels of parallelism allows access to ground and excited state dynamics for large multichromophoric assemblies. The parallel excitonic framework is in good agreement with much more computationally demanding TDDFT calculations of the full assembly.

Model order reduction for a family of linear systems with applications in parametric and uncertain systems

We present an approach which creates a linear reduced order model whose transfer function interpolates at certain given points and approximates at others. We describe how to create the corresponding state space system and explain how this method is used for the simulation of large scale parametric systems or several realizations of an uncertain system. In that case, the interpolation points are the H 2 optimal interpolation points of one particular system. The behavior of the method is illustrated in some numerical examples.

The aggregation and diffusion of asphaltenes studied by GPU-accelerated dissipative particle dynamics

The heavy crude oil consists of thousands of compounds and much of them have large molecular weights and complex structures. Studying the aggregation and diffusion behavior of asphaltenes can facilitate the understanding of the heavy crude oil. In previous studies, the fused aromatic rings were treated as rigid bodies so that dissipative particle dynamics (DPD) integrated with the quaternion method can be used to study asphaltene systems. In this work, DPD integrated with the quaternion method is implemented on graphics processing units (GPUs). Compared with the serial program, tens of times speedup can be achieved when simulations performed on a single GPU. Using multiple GPUs can provide faster computation speed and more storage space for simulations of significant large systems. By using large systems, simulations of the asphaltene–toluene system at extremely dilute concentrations can be performed. The determined diffusion coefficients of asphaltenes are similar to that in experimental studies. At last, the aggregation behavior of asphaltenes in heptane was investigated, and the simulation results agreed with the modified Yen model. Monomers, nanoaggregates and clusters were observed from the simulations at different concentrations.

A parallel implementation of the Stokesian Dynamics method applied to the simulation of colloidal suspensions

We present a strategy for the numerical simulation of polydisperse colloidal suspensions on supercomputers using distributed memory. In order to simulate large colloidal systems we have parallelized the Stokesian Dynamics method and combined it with DLVO theory. Through an efficient parallelization using the message passing system MPI we are able to reduce the computational costs of the problem from O(N2) to O(Pnp,max2), where N is the number of particles in the system, P the number of used processes and np,max=⌈NP⌉, respectively. This allows the simulation of large colloidal systems, as the computational time now not only depends on the number of particles, but also on the number of cores, through which the CPU time can be reduced dramatically.

On The Nature of the Halogen Bond.

The wide-ranging applications of the halogen bond (X-bond), notably in self-assembling materials and medicinal chemistry, have placed this weak intermolecular interaction in a center of great deal of attention. There is a need to elucidate the physical nature of the halogen bond for better understanding of its similarity and differences vis-à-vis other weak intermolecular interactions, for example, hydrogen bond, as well as for developing improved force-fields to simulate nano- and biomaterials involving X-bonds. This understanding is the focus of the present study that combines the insights of a bottom-up approach based on ab initio valence bond (VB) theory and the block-localized wave function (BLW) theory that uses monomers to reconstruct the wave function of a complex. To this end and with an aim of unification, we studied the nature of X-bonds in 55 complexes using the combination of VB and BLW theories. Our conclusion is clear-cut; most of the X-bonds are held by charge transfer interactions (i.e., intermolecular hyperconjugation) as envisioned more than 60 years ago by Mulliken. This is consistent with the experimental and computational findings that X-bonds are more directional than H-bonds. Furthermore, the good linear correlation between charge transfer energies and total interaction energies partially accounts for the success of simple force fields in the simulation of large systems involving X-bonds.

The SpineML toolchain: enabling computational neuroscience through flexible tools for creating, sharing, and simulating neural models

Within computational neuroscience the convergence of large scale computation and biology offers the potential to offer new insight into biological function. Conversely, the simulation of large scale biological systems has the potential to revolutionise computing systems and software which are increasingly showing a trend towards lower power and increased parallelism. In order to achieve large scale simulation, computational neuroscience projects must increase in both scale and complexity. With these increases the problems of managing computational models, sharing these models, and simulating the models rapidly and efficiently become increasingly important to the rate of scientific progress. We present a flexible toolchain that seeks to tackle these issues by providing a standardised pipeline for the specification and simulation of large scale neural models (consisting of rate-coded or point neurons). As the basis of this toolchain we have developed SpineML [1], an XML standard for model description. SpineML ensures compatibility between different models as it contains no simulator or authoring tool specific implementation details, allowing models developed using different tools with different simulators to be combined without alteration. This allows collaborators to develop models independently, using SpineML as a common interchange format. Models authored in SpineML can easily be used in a ‘systems integration’ approach to computational neuroscience. Numerous models can be combined as system components to form models of increasing complexity. This approach can be used to accelerate the progress of research by avoiding the recreation of existing models. SpineML also facilitates model sharing through layer oriented design [2] which provides clear separation of experimental details. This allows changes to simulation parameters, inputs and outputs, and model properties to be performed in an experiment file without changing the model itself and allows collaborators to test different conditions without making changes to the shared model. SpineML provides syntax governed by clearly specified XML schemas. This ensures compatibility between models, however does not place any restriction is on the tools used for authoring models. We have produced a reference model creation tool, SpineCreator, which provides a graphical interface to the SpineML format. This tool allows SpineML models to be created using graphical methods, and visualized in 3D, and is designed to require no knowledge of programming to use, thus making the toolchain accessible to the widest range of users possible. Simulators can support SpineML network descriptions at two levels, with the lowest level interface allowing greater flexibility of the specification and connectivity of the neuronal dynamics. Multi-level support allows SpineML models to be supported and deployed on a wide range of simulators and simulator hardware. Simulation support is currently available using a range of traditional simulators (via PyNN), a reference simulator (BRAHMS) and simulation support for GPUs via GeNN is under development. In addition to this a toolchain for the low power, highly distributed, SpiNNaker hardware architecture has been developed. SpineML facilitates the accessibility of this neuromorphic platform by providing a modeling framework which can be used to develop and test models using traditional simulators which can then be deployed without change directly to neuromorphic hardware. This work was funded by the EPSRC (BIMPA Project (EP/G015740/1) and Green Brain Project (EP/J019690/1)).

Recent advances toward a general purpose linear-scaling quantum force field.

ConspectusThere is need in the molecular simulation community to develop new quantum mechanical (QM) methods that can be routinely applied to the simulation of large molecular systems in complex, heterogeneous condensed phase environments. Although conventional methods, such as the hybrid quantum mechanical/molecular mechanical (QM/MM) method, are adequate for many problems, there remain other applications that demand a fully quantum mechanical approach. QM methods are generally required in applications that involve changes in electronic structure, such as when chemical bond formation or cleavage occurs, when molecules respond to one another through polarization or charge transfer, or when matter interacts with electromagnetic fields. A full QM treatment, rather than QM/MM, is necessary when these features present themselves over a wide spatial range that, in some cases, may span the entire system. Specific examples include the study of catalytic events that involve delocalized changes in chemical bonds, charge transfer, or extensive polarization of the macromolecular environment; drug discovery applications, where the wide range of nonstandard residues and protonation states are challenging to model with purely empirical MM force fields; and the interpretation of spectroscopic observables. Unfortunately, the enormous computational cost of conventional QM methods limit their practical application to small systems. Linear-scaling electronic structure methods (LSQMs) make possible the calculation of large systems but are still too computationally intensive to be applied with the degree of configurational sampling often required to make meaningful comparison with experiment.In this work, we present advances in the development of a quantum mechanical force field (QMFF) suitable for application to biological macromolecules and condensed phase simulations. QMFFs leverage the benefits provided by the LSQM and QM/MM approaches to produce a fully QM method that is able to simultaneously achieve very high accuracy and efficiency. The efficiency of the QMFF is made possible by partitioning the system into fragments and self-consistently solving for the fragment-localized molecular orbitals in the presence of the other fragment’s electron densities. Unlike a LSQM, the QMFF introduces empirical parameters that are tuned to obtain very accurate intermolecular forces. The speed and accuracy of our QMFF is demonstrated through a series of examples ranging from small molecule clusters to condensed phase simulation, and applications to drug docking and protein–protein interactions. In these examples, comparisons are made to conventional molecular mechanical models, semiempirical methods, ab initio Hamiltonians, and a hybrid QM/MM method. The comparisons demonstrate the superior accuracy of our QMFF relative to the other models; nonetheless, we stress that the overarching role of QMFFs is not to supplant these established computational methods for problems where their use is appropriate. The role of QMFFs within the toolbox of multiscale modeling methods is to extend the range of applications to include problems that demand a fully quantum mechanical treatment of a large system with extensive configurational sampling.

Latency-Based Approach to the Simulation of Large Power Electronics Systems

We define a new hybrid simulation scheme that combines latency insertion with nodal and state-space methods. The conditions for stability and multirate time stepping of the hybrid method are defined. Further, we show how to define appropriate flow variable injections so that the hybrid method can be extended to applications involving either multiple technical domains or hardware in the simulation loop. The hybrid method permits significant increases of simulation speed while contributing less than 1% to RMS simulation error.

Quantum mechanical fragment methods based on partitioning atoms or partitioning coordinates.

Conspectus The development of more efficient and more accurate ways to represent reactive potential energy surfaces is a requirement for extending the simulation of large systems to more complex systems, longer-time dynamical processes, and more complete statistical mechanical sampling. One way to treat large systems is by direct dynamics fragment methods. Another way is by fitting system-specific analytic potential energy functions with methods adapted to large systems. Here we consider both approaches. First we consider three fragment methods that allow a given monomer to appear in more than one fragment. The first two approaches are the electrostatically embedded many-body (EE-MB) expansion and the electrostatically embedded many-body expansion of the correlation energy (EE-MB-CE), which we have shown to yield quite accurate results even when one restricts the calculations to include only electrostatically embedded dimers. The third fragment method is the electrostatically embedded molecular tailoring approach (EE-MTA), which is more flexible than EE-MB and EE-MB-CE. We show that electrostatic embedding greatly improves the accuracy of these approaches compared with the original unembedded approaches. Quantum mechanical fragment methods share with combined quantum mechanical/molecular mechanical (QM/MM) methods the need to treat a quantum mechanical fragment in the presence of the rest of the system, which is especially challenging for those parts of the rest of the system that are close to the boundary of the quantum mechanical fragment. This is a delicate matter even for fragments that are not covalently bonded to the rest of the system, but it becomes even more difficult when the boundary of the quantum mechanical fragment cuts a bond. We have developed a suite of methods for more realistically treating interactions across such boundaries. These methods include redistributing and balancing the external partial atomic charges and the use of tuned fluorine atoms for capping dangling bonds, and we have shown that they can greatly improve the accuracy. Finally we present a new approach that goes beyond QM/MM by combining the convenience of molecular mechanics with the accuracy of fitting a potential function to electronic structure calculations on a specific system. To make the latter practical for systems with a large number of degrees of freedom, we developed a method to interpolate between local internal-coordinate fits to the potential energy. A key issue for the application to large systems is that rather than assigning the atoms or monomers to fragments, we assign the internal coordinates to reaction, secondary, and tertiary sets. Thus, we make a partition in coordinate space rather than atom space. Fits to the local dependence of the potential energy on tertiary coordinates are arrayed along a preselected reaction coordinate at a sequence of geometries called anchor points; the potential energy function is called an anchor points reactive potential. Electrostatically embedded fragment methods and the anchor points reactive potential, because they are based on treating an entire system by quantum mechanical electronic structure methods but are affordable for large and complex systems, have the potential to open new areas for accurate simulations where combined QM/MM methods are inadequate.

Voxel based parallel post processor for void nucleation and growth analysis of atomistic simulations of material fracture

Molecular dynamics (MD) simulations are used in the study of void nucleation and growth in crystals that are subjected to tensile deformation. These simulations are run for typically several hundred thousand time steps depending on the problem. We output the atom positions at a required frequency for post processing to determine the void nucleation, growth and coalescence due to tensile deformation. The simulation volume is broken up into voxels of size equal to the unit cell size of crystal. In this paper, we present the algorithm to identify the empty unit cells (voids), their connections (void size) and dynamic changes (growth and coalescence of voids) for MD simulations of large atomic systems (multi-million atoms). We discuss the parallel algorithms that were implemented and discuss their relative applicability in terms of their speedup and scalability. We also present the results on scalability of our algorithm when it is incorporated into MD software LAMMPS.

Hydrodynamic shock wave studies within a kinetic Monte Carlo approach

We introduce a massively parallelized test-particle based kinetic Monte Carlo code that is capable of modeling the phase space evolution of an arbitrarily sized system that is free to move in and out of the continuum limit. Our code combines advantages of the DSMC and the Point of Closest Approach techniques for solving the collision integral. With that, it achieves high spatial accuracy in simulations of large particle systems while maintaining computational feasibility. Using particle mean free paths which are small with respect to the characteristic length scale of the simulated system, we reproduce hydrodynamic behavior. To demonstrate that our code can retrieve continuum solutions, we perform a test-suite of classic hydrodynamic shock problems consisting of the Sod, the Noh, and the Sedov tests. We find that the results of our simulations which apply millions of test-particles match the analytic solutions well. In addition, we take advantage of the ability of kinetic codes to describe matter out of the continuum regime when applying large particle mean free paths. With that, we study and compare the evolution of shock waves in the hydrodynamic limit and in a regime which is not reachable by hydrodynamic codes.

Shock-induced plasticity and the Hugoniot elastic limit in copper nano films and rods

Shock deformation of copper nano-films and nano-rods is examined with Molecular Dynamics (MD) simulations. The influence of the small system size on the onset of plasticity, its origin resulting from the nucleation of dislocation loops, and its reversible nature are determined. While simulations of large systems with periodic boundary conditions indicate that tremendous axial stresses are needed to induce plastic deformation in perfect copper crystals, the present results suggest that the stress levels needed to initiate irreversible plasticity in nano-rods are more than one order of magnitude smaller than what has been reported for bulk single crystals. MD studies of nano-films show that shock waves are purely elastic up until the Hugoniot elastic limit of PHEL ≈ 30–40 GPa, at which point Shockley partial dislocations are internally nucleated at the shock front. However, our recent experiments on shocked nano-rods show that plasticity is evident at much lower axial stress levels, on the order of 1–2 GPa. The present MD simulations of shocked nano-rods show that Shockley partial dislocations prefer to nucleate at lower stresses from the rod surface, at PHEL ≈ 1–2 GPa, consistent with our concurrent experimental observations, leading to surface step formation and mechanical damage. Nucleated dislocations are found to be Shockley partials in the [100] and [111] oriented nano-rods, with the additional presence of perfect dislocations in the latter. MD simulations of rarefaction shock waves in nano-films indicate that they can be spalled via a mechanism of nano-void nucleation, growth and coalescence at the spall plane. The origin of these nano-voids is shown to be at the intersection of stacking faults on conjugate slip {111} planes. Spallation by void nucleation and coalescence is found not to be achievable in nano-rods. Rarefaction shocks with high stresses were found to either severely deform or melt the nano-rod before it can be spalled.

A New Strategy for Coarse-Grained Protein Simulations: Smoothed Energy Tables

Coarse grained simulation methods allow the simulation of large biomolecular systems with less computational expense than the corresponding all-atom simulations. The savings come from two sources: faster computation due to a reduced number of particles, and improved sampling due to a smoother free energy surface. Many commonly used coarse-grained models suffer from serious limitations, such as being unable to properly model protein secondary structure without the addition of unphysical restraints. We have constructed a novel coarse-grained Monte Carlo method based on dividing proteins into nearly rigid fragments, constructing distance and orientation-dependent tables of the interaction energies between those fragments, and applying potential energy smoothing techniques to those tables. Preliminary results on peptides indicate that the new method is able to preserve α-helices without additional restraints. In addition, when sufficient smoothing is applied, the new method also shows an improvement in sampling per unit computation time compared to Monte Carlo simulation without tabulation or atomistic molecular dynamics.

Frequency Weighted Model Order Reduction Technique and Error Bounds for Discrete Time Systems

Model reduction is a process of approximating higher order original models by comparatively lower order models with reasonable accuracy in order to provide ease in design, modeling and simulation for large complex systems. Generally, model reduction techniques approximate the higher order systems for whole frequency range. However, certain applications (like controller reduction) require frequency weighted approximation, which introduce the concept of using frequency weights in model reduction techniques. Limitations of some existing frequency weighted model reduction techniques include lack of stability of reduced order models (for two sided weighting case) and frequency response error bounds. A new frequency weighted technique for balanced model reduction for discrete time systems is proposed. The proposed technique guarantees stable reduced order models even for the case when two sided weightings are present. Efficient technique for frequency weighted Gramians is also proposed. Results are compared with other existing frequency weighted model reduction techniques for discrete time systems. Moreover, the proposed technique yields frequency response error bounds.