\section{\label{sec:level1}Introduction}
\indent Neural network model chemistries (NNMCs) greatly reduce the computational effort needed to simulate chemical systems with ab-initio accuracy\cite{Snyder:2013aa,brockherde2017bypassing,snyder2012finding,li2016understanding,li2016pure,vu2015understanding,handley2010potential,chmiela2017machine,behler2011neural,behler2007generalized,shakouri2017accurate,behler2017first, han2017deep,yao2017many,yao2016kinetic,yao2017intrinsic,khaliullin2011nucleation,bartok2010gaussian, mones2016exploration,gastegger2017machine,kobayashi2017neural, carpenter2017empirical,kolb2017discovering,kruglov2017energy, lubbers2017hierarchical,mills2017deep,wu2017internal,khorshidi2016amp,shao2016communication,zhang2014effects, li2015permutationally2,medders2015representation,medders2013critical,reddy2016accuracy,riera2017toward,moberg2017molecular, conte2015permutationally,manzhos2014neural, manzhos2009fitting,malshe2010input,peterson2016acceleration, khorshidi2016amp,piquemal2017preface,cubuk2017representations, MBCoarseCsanyi,FFparafit,li2017machine}. They can be used to predict molecular properties \cite{CM, BoBs,lopez2014modeling, pilania2013accelerating, schutt2014represent, ma2015machine, nelson2012nonadiabatic, janet2017predicting, janet2017resolving, yao2017intrinsic, hase2017machine, mcgibbon2017improving, bereau2017non, grisafi2017symmetry, isayev2017universal, ghiringhelli2017learning, ouyang2017sisso, faber2017prediction, schutt2017quantum, MSLSTM, timoshenko2017supervised}, and design new materials or drugs\cite{li2017learning, ramsundar2017multitask, hachmann2011harvard,hachmann2014lead, isayev2015materials, kim2017machine, segler2017towards,olivares2011accelerated,guimaraes2017objective,wei2016neural, gomez2016automatic,jinnouchi2017predicting,ulissi2017machine,sun2017machine}. In spite of growing popularity, most neural network methods are still only used by their developers and are customized for a single application. The paucity of literature describing transferable accuracy, confusion about what physics can be reproduced, and dearth of open software is slowing adoption. This paper develops an open-source, transferable neural network model chemistry. We show that NNMCs are easily hybridized with physical contributions to molecular energies familiar from Molecular Mechanics and corrections to Density Functional Theory(DFT)\cite{grimme2006semiempirical}. This approach combines the best of both worlds, yielding predictable reproduction of physical long-range forces, but also featuring a linear-scaling inductive charge model which is cheaper than a Thole model\cite{thole1981molecular} yet more accurate than fixed charges. \\ 
\indent     Our group is one of several who have been pursuing transferable and black-box neural network model chemistries\cite{yao2017intrinsic, brockherde2017bypassing, smith2017ani, medders2013critical, behler2007generalized, han2017deep}. The field is growing so rapidly, that most non-practitioners cannot keep up with the capabilities of existing models and the outstanding problems. Readers may not appreciate that a model can achieve chemical accuracy for energies but have uselessly noisy forces. Models which provide energies at equilibrium, and those treating a fixed molecule or stoichiometry are now reliably produced\cite{han2017deep}. We will show that TensorMol-0.1 yields usefully accurate predictions of forces out-of-equilibrium by showing reproduction of infrared spectra which closely approximate our source model chemistry (wB97X-D, 6-311G**)\cite{chai2008long}, and molecular dynamics. We outline several tricks which are required to ensure the stability of long-time molecular dynamics.\\
\indent     Another distinguishing feature of our approach is the open-source nature of our package, which uses the TensorFlow tensor algebra system to compute descriptors and forces. Rather than a monolithic black-box TensorMol-0.1 is a modular collection of differentiable chemical models written in concise Python code. The components are easily joined together and extended. The methodology can be used to propagate dynamics for large molecules ($10^5$ atoms) with GPU acceleration on simple laptop computers. No significant expertise, force field refinement, or other interventions are needed to apply the method to a molecule of interest, so long as the elements are supported. The package is also interfaced with the I-PI path integral package\cite{ceriotti2014pi}, to allow for quantum simulations and enhanced sampling. \\

\section{\label{sec:level1}Methods}

\begin{figure}[t]
\includegraphics[width=0.45\textwidth]{scheme.pdf}
\caption{The schematic graph of TensorMol-0.1. Each element has its own charge network and energy network. The charge network predicts atomic charges that yield the ab-initio dipole moment. The Behler-Parinello type energy network produces a short-range embedded atomic energy, which is summed with the electrostatic energy and Van-Der Waals energy to predict the total atomization energy of molecules at and away from equilibrium. The whole computation is included in the TensorFlow graph. Analytical GPU-accelerated forces are symbolically generated by a single line of code.}
\label{Fig1}
\end{figure}

\indent     The community of Neural-Network model chemistry developers is rapidly improving the accuracy and generality of these reactive force fields\cite{behler2017first,  shakouri2017accurate, bartok2010gaussian, deringer2017machine, li2015permutationally2, smith2017ani}. The model of this paper includes several components which were the subject of recent developments in other groups\cite{behler2007generalized,behler2017first,morawietz2013density, smith2017ani, gastegger2017machine}. We will describe the details here from the bottom up citing prior art. Our notational convention will be that $i,j,k...$ are indices of atoms, $q_i$ is the charge on atom i, $z,x,y$ are atomic numbers, $\mathcal{A},\mathcal{B}, \mathcal{C}$ are molecules, and $\alpha, \beta... $ are indices of basis functions which are a product of radial and angular functions. If a function depends on all the atomic coordinates of a molecule it will be written as a vector, and those which depend on only few will be given explicit indices. The energy of TensorMol-0.1 is expressed as a sum of a short-range embedded N-body potential\cite{behler2017first}, and long-range electrostatic potential and Van-Der-Waals force: 
\begin{eqnarray}
E(\vec{R}) =&&\sum_i E^\text{BP}_{z_i} (S_\alpha(\vec{R})) +\sum_{ij} E^\text{DSF}_{ij}(q_{ij}(S_\alpha(\vec{R})),R_i,R_j)\nonumber\\&&+E^\text{VDW}(\vec{R}_{ij})
\end{eqnarray}
In the above expression $E_{z_i}$ is a Behler-Parinello type energy network for the element $z$ for atom $i$. This n-body potential takes as its argument $S_\alpha$ the modified symmetry functions of Isayev and coworkers\cite{smith2017ani}: 
\begin{eqnarray}
S_\alpha (radial) = && \sum_{j\neq i}e^{{-\eta(R_{ij}-R_{s})^2}}f_c(R_{ij}) \\
S_\alpha (angular) =&& 2^{1-\zeta }\sum_{j\neq i, j\neq k}(1+cos(\theta_{ijk}-\theta_s))^{\zeta}\nonumber\\&&\times e^{-\eta(\frac{R_{ij}+R_{ik}}{2}-R_s)^2}f_c(R_{ij})f_c(R_{ik})
\end{eqnarray}
Modern machine learning frameworks provide automatic differentiation of tensor algebraic expressions, allowing a force-field developer to obtain the gradient of a molecular potential $\frac{dE(\vec{R})}{d\vec{R}}$ in a single line of code, once the expression for $E(\vec{R})$ has been written. An important feature of our code is that this symmetry function is coded within the TensorFlow system\cite{tensorflow2015-whitepaper}, so all the parameters of this descriptor can be variationally optimized alongside the network weights. Our implementation of the symmetry function employs a list of nearest-pairs and triples within radial cutoffs such that the scaling of the overall network is asymptotically linear. On an ordinary laptop equipped with only a CPU a force/energy call on 20,000 atoms is less than a minute. \\
\indent     The second term of our energy expression is the Damped-Shifted Force (DSF) Coulomb energy of Gezelter and coworkers\cite{fennell2006ewald}. The charges are obtained from a sub-network which reproduces molecular dipole moments. Our charge model enforces conservation of total charge by evenly spreading any required neutralizing charge over the entire molecule or unit cell. The Damped-shifted force ensures long range continuity and differentiability of the effective Coulomb potential with smooth cutoffs. We modify the DSF kernel at short range with an "elu" type non-linearity, such that the forces within the radius of the Behler-Parinello symmetry function smoothly approach zero avoiding singularities and interference with the Behler-Parinello many-body potential.
\begin{eqnarray}
E_{DSF}= &&\left\{\begin{matrix}  E_\text{DSF(Original)}  \hfill   R > R_\text{switch}
\\ q_iq_j(a_\text{elu}e^{R-R_\text{switch}}+\beta_\text{elu})\hfill R< R_\text{switch}
\end{matrix}\right.
\end{eqnarray}
where $E_\text{DSF(Original)}$ is the energy of DSF kernel\cite{fennell2006ewald},  $R_{switch}$ is the short range cutoff for the "elu" kernel. $\alpha_\text{elu}$ and $\beta_\text{elu}$ are chosen so that the value and the gradient of $E_\text{DSF}$ are continuous at $R_\text{switch}$.
The damped-shifted force is well-suited to combination with Neural Network models because it requires no Fourier transformation to treat periodic systems with linear scaling and maps well onto TensorFlow. The last term is the van-der waals energy, which is calculated by following Grimme's C6 scheme \cite{grimme2006semiempirical}.\\
\indent We employed a two step training approach. First, the charge networks are trained to learn the atom charges that predict the dipole moment. The loss function can be written as following:
\begin{align}
 L_{dipole} = &&\sum_\mathcal{A} (\frac{\mu^\text{DFT}_\mathcal{A} - \mu^{NN}_\mathcal{A}(q_i, q_j,...)}{N_\text{atom}})^2
\end{align}
After the charge training is converged, we train the energy network. During the energy network training, the weights in charge networks are kept frozen, but they are still evaluated to calculate the electrostatic energy that is added to construct the total energy. Our Behler-Parinello many-body potential also absorbs the shape of the transition between many-body and electrostatic regions. The learning target for the energy network includes both the DFT energy and DFT force. The loss function for the energy network training is:
\begin{align}
 L_\text{energy} =&& \sum_\mathcal{A} (\frac{E^\text{DFT}_\mathcal{A} - E^\text{NN}_\mathcal{A}}{N_{atom}})^2 + \gamma \sum_\mathcal{A} (\frac{F^\text{DFT}_\mathcal{A} - F^\text{NN}_\mathcal{A}}{N_\text{atom}})^2
\end{align}
where $E^\text{NN}$ is obtained according to equation 1, $F^\text{NN}$ is calculated by taking the gradient of $E^\text{NN}$ respect to the coordinates of the atoms. $N_\text{atom}$ is the number of the atoms in the system and $\gamma$ is a parameter that controlling the portion of force loss. We employ $\gamma=0.05$. 
\indent  We trained two neural networks based on two sets of data. One network ("water network") was trained on a database that includes $\sim$370,000 water clusters with 1 water molecule to 21 water molecules. The other network was trained on $\sim$3,000,000 different geometries of $\sim$ 15,000 different molecules that only contains C, H, O and N and up to 35 atoms. Since these 15K molecules were sampled randomly from the chemspider database, therefore we will refer this network as "chemspider network" in the following text. The training data were sampled using metadynamics and calculated by Qchem package\cite{shao2015advances} with WB97X-D\cite{chai2008long},  exchange correlation functional and 6-311G** basis set.\\

Each charge network and energy network contains three fully-connected hidden layers with 500 hidden neurons in each layer. For chemspider network, a network with three hidden layers with 2000 hidden neurons in each layers is used for each charge network and energy network. L2 regularization and dropout \cite{srivastava2014dropout} on last layer were used to prevent overfitting with a dropout probability of 0.3. We chose a softplus function as the non-linear activation function after extensive experimentation, and used the Adaptive moment quasi-Newton solver (Adam)\cite{kingma2014adam} to fix the weights of the network. The test sets were separated from training data by choosing a random 20\% of molecules at the outset which were kept independent throughout. Besides water we will present calculations from molecules strictly absent from either the training or test set. 


\indent     To obtain scalable efficiency TensorMol uses neighborlists within cutoffs. This allows double precision energy, charge, force calculations of up to 24,000 atoms to execute in less than 90 seconds on a 2015 Intel i7 2.5GHz MacBook pro (Fig. \ref{fig:timings}). Periodic evaluations are achieved by tessellation of a unit cell with summation of energies for atoms within the cell. This results in  roughly a factor of three overhead in computational time. Speedups are obtained automatically for computers with GPUs or single-precision calculations. 
\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{timing_water.pdf}
    \caption{Aperiodic timings of an energy, charge, force call for cubic water clusters at a density of 1 gm/cm$^3$. The largest $\sim60$ Angstrom cube is 4x larger than the electrostatic cutoff. The slope of a log-log version of this curve is near unity, indicating the wall-time scaling of TensorMol.}
    \label{fig:timings}
\end{figure}


\section{\label{sec:level1}Results}
The root mean square error (RMSE) on the independent test set of the energy is 0.054 kcal/mol per atom and the RMSE of the force is 0.49 kcal/mol/\AA. The left panel of figure \ref{fig:trimer} plots the potential energy surface (PES) of a water trimer when one of the water is pulled away from the other two. One can see our neural network PES is not only in good agreement with the PES of target method but also smooth. To achieve this we use a variation of the soft-plus neuron rather than the rectified linear units which are popular in computer science. The latter train more efficiently, but produce discontinuous forces. \\
\indent     The right panel shows the fraction of each of the three energy components in equation 1 to the binding energy along the trimer dissociation coordinate.  At short range, most of the binding energy is contributed by the N-body neural network potential. When the distance between the monomer and the dimer approach the cutoff distance of the neural network, the contribution of neural network potential starts to decrease and the contribution of electrostatic potential increases. After 6\ \AA\ where the neural network symmetry functions on the atoms in the monomer have no contribution from the dimer, the neural network force drops smoothly to zero and the electrostatic interaction dominates. The small difference in the energy at 7\ \AA\ is due to the difference between the Madelung energy given by the learned charges, and the genuine physical cohesive force at this distance. The dimer and monomer are beyond the symmetry function sensory radius, and so the charges are constant in this region. Future iterations of the charge network will use local-field information to improve this region of the PES. The learned inductive charges are of high quality considering their linear scaling cost. Figure \ref{fig:dimer} shows the PES and dipole change of a water dimer when the hydrogen bond is broken by rotating the OH bond. Both the PES and dipole change fit well with the DFT results. \\

\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{trimerbinding.pdf}
    \caption{Left panel: PES of water trimer when one water is pulled away from the other two. Right panel: Percentage contribution of binding energy between the water that is pulled away and the other two water from Behler-Parrinello atom-wise energy, electrostatic energy and van-der waals energy. Behler-Parrinello atom-wise energy contribute to most of the binding energy at the short range and electrostatic energy is the dominant contribution at long range.}
    \label{fig:trimer}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{dimerdipole.pdf}
    \caption{Top left panel: PES of breaking a hydrogen bond between two water by rotating one water around O-H bond. Top right, bottom left and bottom right panels: change of x, y, z, dipole component during the rotation, respectively}
    \label{fig:dimer}
\end{figure}

Given the increased dimension of the Hessian, it is naturally a more stringent test to reproduce forces and infrared spectra than it is to simply produce energies. The left panel and right panel of figure \ref{fig:waterir} show the optimized geometries and IR spectra of a 10 water cluster and 20 water cluster generated with our force field and DFT, respectively. Each method uses its own equilibrium geometry, so this also tests TensorMol-0.1's reproduction of non-covalent geometry. The RMSE of the distance matrix between DFT optimized geometry and TensorMol optimized geometry are 0.062\ \AA\ for the 10 water cluster and 0.180\ \AA\ for the 20 water cluster. Our force field quantitatively reproduces the DFT IR both in terms of frequencies and intensities, especially for the water bend modes and inter-monomer modes. The Mean Absolute Error (MAE) of frequencies at that those two regions are 33.2 cm\textsuperscript{-1} for the 10 water cluster and 16.2 cm\textsuperscript{-1} for the 20 water cluster. The error is slightly larger at  water OH stretching region with a MAE of 34.2 cm\textsuperscript{-1} and 13.1 cm\textsuperscript{-1}, respectively. This accuracy is comparable to high quality polarizable water force fields\cite{medders2015representation}.\\

\begin{table}
\caption{\label{tab:table1} Training details and test RMSE of each learning target. The unit of energy RMSE, gradient RMSE and dipole RMSE is kcal/mol per atom, kcal/mol/\AA\ per atom and Debye per atom, respectively.}
\begin{ruledtabular}
\begin{tabular}{lcc}
&Water Network&Chemsipider Network\\
\hline
Num of training case & 370844 & 2979162\\
Training time (days)\footnote{Training was done on single Nvidia K40 GPU}& 3  & 10 \\
Energy RMSE& 0.054 & 0.24\\
Gradient RMSE& 0.49 & 2.4\\
Dipole RMSE& 0.0082  & 0.024\\
\end{tabular}
\end{ruledtabular}
\end{table}

\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{water_IR.pdf}
    \caption{Simulated harmonic IR spectrum of 10 water cluster (left panel) and 20 water cluster (right panel) generated by WB97X-D/6-311G**(dashed orange line) and TensorMol force field (solid blue line).}
    \label{fig:waterir}
\end{figure}
Compared with traditional force fields, one major advantage of TensorMol is its reactivity. TensorMol is able to simulate a concerted proton transfer in a water hexamer, finding a minimum energy transition path. The PES's calculated by nudged elastic band (NEB) method \cite{henkelman2000climbing} with the TensorMol force field and DFT are shown in figure \ref{fig:waterneb}. The barrier height predicted by TensorMol is 36.3 kcal/mol, which is 6.7 kcal/mol lower than the prediction of DFT, which is remarkable considering the dearth of transition structures in the training data. Our sampling of molecular geometries uses a meta-dynamics procedure described elsewhere, so these proton transfers occur in the training data although extremely infrequently. 

Encouraged by our water results, we developed a force field with applicability across the chemical space spanned by CNOH. The Chemspider dataset that we used to train our force field covers a vast chemical space containing 15 thousand different molecules and 3 millions geometries. The geometries are generated using a meta-dynamics procedure\cite{barducci2008well}, which ensures that each new geometry is a fresh part of chemical space, energies up to 400$k_bT$ are sampled in the data. We describe the details of this meta-dynamics sampling algorithm, which we have found vital to achieving robust and transferrable force-fields elsewhere. The diversity of structures makes learning the chemspider dataset a much harder task for neural networks, the test set RMSE of energy is 0.24 kcal/mol per atom and RMSE of force is 2.4 kcal/mol per atom. More importantly, the model usefully reproduces several elements of molecular structure at and away from equilibrium for molecules outside its training set. It robustly optimizes the geometries of typical organic molecules to structures that match DFT well, and yields infrared frequencies and intensities in good agreement with ab-initio calculations. It is a black-box method which does not rely on any atom type, connectivity, etc as one would need to specify in a traditional classical force-field. The few proteins we have examined remain stable and near their experimental structures when optimized or propagated at room temperature using the TensorMol-0.1 force field.\\
\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{water_neb.pdf}
    \caption{Reaction energy profile converged from a Nudged elastic band along the reaction coordinate of conservative proton transfer in a water hexamer cluster.}
    \label{fig:waterneb}
\end{figure}
\indent Morphine is not included in our training set. The left panel of figure \ref{fig:morphine} shows the geometry of morphine that is optimized with our force field. The RMSE of bond lengths predicted by our forcefield is 0.0067\ \AA\ and the RMSE of angles is 1.04 degrees compared with the source DFT model chemistry. The right panel plots the harmonic IR spectra generated by each method at their respective optimized geometries. One can see IR spectrum from our force field is in good agreement with the DFT IR spectrum. The MAE of our force field frequencies is 13.7 cm\textsuperscript{-1} compared with DFT frequencies. Figure \ref{fig:chemspider} shows comparisons of IR spectrum that are generated from these two methods for aspirin, typrosine, caffeine and cholesterol. All these four molecules are not included in the training set. The MAE of the frequencies predicted by our field is less than 20 cm\textsuperscript{-1} for all the four molecules compared with target DFT frequencies. The concept of a chemical bond and force constant are not enforced in any way, yet good agreement with DFT is obtained at a tiny fraction of the original cost. \\
\indent    Traditional harmonic vibrational spectra require quadratic computational effort, which works against the speed advantage of a NNMC. For large systems one can use the molecular dynamics functionality of TensorMol to simulate infrared spectra, Fourier transforming the dipole-dipole correlation function of conservative Newtonian dynamics whose cost grows linearly with the size of the system. The lower left panel of Figure \ref{fig:morphine} shows the same infrared spectrum produced by propagation in TensorMol-0.1, also showcasing the good energy conservation of TensorMol. Unlike a traditional force-field it's non-trivial to obtain smoothly differentiable NNMC's. 64-bit precision needs to be used the network cannot be made too flexible and smooth versions of typical rectified linear units need to be used. Our package can be used in this way to simulate IR of large systems with linear cost.\\
\indent     TensorMol-0.1 has a relatively simple treatment of electrostatic and Van-Der-Waals forces which we would like to augment in the future with a many-body dispersion scheme\cite{tkatchenko2012accurate}. However a main adantage of TensorMol-0.1's approach is its very low cost. No self-consistent polarization equation is solved even though the charges are inductive, and so it is easy to inexpensively calculate the electrostatic energies of even very large molecules. At shorter ranges, non-covalent interactions like hydrogen bonds are dealt with by the Behler-Parinello portion of the network. The Chemspider training data includes some examples of dimers and intra-molecular hydrogen bonds. To our surprise the treatment of inter-molecular interactions which were not targets for TensorMol-0.1 are satisfactory. Figure \ref{fig:DNA} shows the optimized geometries and binding energies of two DNA base pairs calculated by our force field. The target DFT method predicts a binding energy of 18.3 kcal/mol for the thymine-adenine (TA) pair and a binding energy of 32.4 kcal/mol for the guanine-cytosine (GC) pair. The prediction of our force field is 1.2 kcal/mol less for TA pair and 2.0 kcal/mol larger for GC pairs relative to DFT. \\ 
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{morphine_IR.pdf}
    \caption{Morphine geometry that is optimized by TensorMol-0.1 (upper right panel) and its harmonic IR spectrum simulated by WB97X-D/6-311G**(dashed orange line) and TensorMol force field (solid blue line) (upper left panel). Lower panels show TensorMol's real-time IR spectrum vs. DFT (left) and the conservation of energy maintained by the smoothness of the energy (right).}
    \label{fig:morphine}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{chemspider_IR.pdf}
    \caption{Harmonic IR spectrum of four different molecules simulated by WB97X-D/6-311G**(dashed orange line) and TensorMol-0.1. All the molecules are not included in the training set.}
    \label{fig:chemspider}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{DNA.pdf}
    \caption{Binding energy between the DNA base pairs vs. $\omega$B97x-D with methods at their optimized geometries. The difference between DFT and TensorMol binding energy is < 2 kcal/mol.}
    \label{fig:DNA}
\end{figure}
\indent     One holy grail of the field of Neural Network model chemistries is to simulate biological chemistry without QM-MM or bespoke force-fields. Protein simulation also demonstrates several important features of a neural network model chemistry: reasonable inter-molecular forces, stability, scalability and generalization far from small-molecule training data. TensorMol-0.1 was not trained on any peptide polymers and includes no biological data of any sort. To our pleasant surprise, even this first iteration of Neural Network model chemistry is accurate enough to perform rudimentary studies of small proteins. Figure \ref{fig:protein} shows geometries sampled from a 1 picosecond, periodic, 300K TensorMol dynamics NVT trajectory in explicit solvent. The initial structure (included in the supplement) was generated from the PDB structure 2MZX using OpenMM's automatic solvation and hydrogenation scripts\cite{OpenMM}, but includes nothing but atom coordinates. This short alpha-helix is stable, both in optimizations and dynamics, and the structures sampled during the dynamics superficially resemble the solution NMR structure. Traditional force fields will always be less expensive (by some prefactor) than NNMCs, yet the reactivity advantages of NNMCs and the ease of set up will probably lead to a rapid adoption of these methods in the biological community. 
\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{protein.pdf}
    \caption{Left panel shows samples from a 1 picosecond NVT (Nos\'e) trajectory of solvated 2MZX at 300K simulated by our TensorMol force field in explicit water. Right panel is the NMR structure of 2MZX from the PDB database.}
    \label{fig:protein}
\end{figure}
\subsection{Discussion and Conclusions}
\indent     We have presented a transferable neural network model chemistry with long-range Coulombic physics, and a short-range n-Body potential, TensorMol-0.1. The model is integrated in a concise open-source python package which provides many of the types of simulation commonly used in chemistry. The method can be used to scan conformational and chemical space along the singlet neutral potential energy surface with high throughput and accuracy using nothing but atomic coordinates. TensorMol-0.1 is not the final iteration of Neural-Network model chemistry, although it shows just how easily DFT-quality predictions can be made by models with drastically lower cost. Inexpensive post-DFT corrections such as Many-Body Dispersion \cite{tkatchenko2012accurate} will become even more powerful when integrated with these potentials, opening the door to quantitative treatments of large systems. These methods will compete aggressively with DFT packages, and provide an interesting complement to QM-MM type simulations in the near future. \\
\indent     Methods like TensorMol-0.1 suggest the following directions for the future improvement of Neural-Network model chemistries:
\begin{itemize}
\item Generalize descriptors to encode other physical atom properties besides charge (spin, polarizability)
\item Develop accurate descriptors whose cost grows linearly with the number of elements treated
\item Extend the range of the N-Body embedding
\item Explore the hierarchy of physical detail between Force-Fields and semi-empirical electronic structure. 
\end{itemize}
These goals must be pursed alongside honest test suites with open programs and data. NNMCs that cannot produce forces and MD trajectories with a demonstrable and compelling cost advantage over DFT should only be pursued if they offer interesting qualitative insights. 

\begin{acknowledgments}
The authors gratefully acknowledge Notre Dame's College of Science for startup funding, Oak Ridge national laboratory for a grant of supercomputer resources and NVidia corporation.
\end{acknowledgments}




\input{tensormol.bbl}

\clearpage

\end{document}
