[bug] documentation not shown for some symbols in latex

[cridemichel]

Issue Prelude

Please complete these steps and check these boxes (by putting an x inside the brackets) before filing your issue:

• [X] I have read and understood YCM's [CONTRIBUTING][cont] document.
• [X] I have read and understood YCM's [CODE_OF_CONDUCT][code] document.
• [X] I have read and understood YCM's [README][readme], especially the [Frequently Asked Questions][faq] section.
• [X] I have searched YCM's issue tracker to find issues similar to the one I'm about to report and couldn't find an answer to my problem. ([Example Google search.][search])
• [X] If filing a bug report, I have included the output of vim --version.
• [X] If filing a bug report, I have included the output of :YcmDebugInfo.
• [X] If filing a bug report, I have attached the contents of the logfiles using the :YcmToggleLogs command.
• [X] If filing a bug report, I have included which OS (including specific OS version) I am using.
• [X] If filing a bug report, I have included a minimal test case that reproduces my issue, using vim -Nu /path/to/YCM/vimrc_ycm_minimal, including what I expected to happen and what actually happened.
• [X] If filing a installation failure report, I have included the entire output of install.py (or cmake/make/ninja) including its invocation
• [X] I understand this is an open-source project staffed by volunteers and that any help I receive is a selfless, heartfelt gift of their free time. I know I am not entitled to anything and will be polite and courteous.
• [X] I understand my issue may be closed if it becomes obvious I didn't actually perform all of these steps.

Thank you for adhering to this process! It ensures your issue is resolved quickly and that neither your nor our time is needlessly wasted.

Issue Details

Dear developers, in latex file if I type "eqnarray" from coc.nvim plugin I get the following documentation in a preview popup window:

however I get no documentation when using YCM, although from other latex symbols I get something like in the following example:

Actually I am using texlab which I have set up as follows in my .vimrc file:

let g:ycm_language_server = [
  \   {
  \     'name': 'latex',
  \     'cmdline': [ 'texlab' ],
  \     'filetypes': [ 'tex' ],
  \     'root_files': [ '.criprojroot' ],
  \   },
  \ ]

to reproduce the problem you can open the following latex file

\documentclass[aps,pre,superscriptaddress]{revtex4-1}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{color}
\begin{document}
\title{Life, the Universe and Everything }
\author{Douglas Adams}
\begin{abstract}
42 is the answer
\end{abstract}
\date{\today}
\maketitle
\section{Introduction}
bla bla
\begin{eqnarray}
  x&=&y\\
  y&=&2
\end{eqnarray}
\end{document}

with the following minimal vimrc file:

let maplocalleader=' '
set nocompatible              " be iMproved, required
call plug#begin('~/.vim/plugged')
Plug 'Valloric/YouCompleteMe' 
call plug#end()            " required
set completeopt=longest,menuone,popup
set completepopup=border:off
if !exists('g:ycm_semantic_triggers')
   let g:ycm_semantic_triggers = {}
endif
let g:ycm_add_preview_to_completeopt=1
let g:ycm_language_server = [
  \   {
  \     'name': 'latex',
  \     'cmdline': [ 'texlab' ],
  \     'filetypes': [ 'tex' ],
  \     'root_files': [ '.criprojroot' ],
  \   },
  \ ]
set noshowmode
let g:ycm_log_level = 'debug'

Diagnostic data

Output of vim --version

VIM - Vi IMproved 8.2 (2019 Dec 12, compilato May 02 2021 18:16:24)
Versione macOS - x86_64
Patch incluse: 1-2825
Compilato da Homebrew
Versione gigante senza GUI.  Funzionalità incluse (+) o escluse (-):
+acl               -farsi             +mouse_sgr         +tag_binary
+arabic            +file_in_path      -mouse_sysmouse    -tag_old_static
+autocmd           +find_in_path      +mouse_urxvt       -tag_any_white
+autochdir         +float             +mouse_xterm       -tcl
-autoservername    +folding           +multi_byte        +termguicolors
-balloon_eval      -footer            +multi_lang        +terminal
+balloon_eval_term +fork()            -mzscheme          +terminfo
-browse            +gettext           +netbeans_intg     +termresponse
++builtin_terms    -hangul_input      +num64             +textobjects
+byte_offset       +iconv             +packages          +textprop
+channel           +insert_expand     +path_extra        +timers
+cindent           +ipv6              +perl              +title
-clientserver      +job               +persistent_undo   -toolbar
+clipboard         +jumplist          +popupwin          +user_commands
+cmdline_compl     +keymap            +postscript        +vartabs
+cmdline_hist      +lambda            +printer           +vertsplit
+cmdline_info      +langmap           +profile           +virtualedit
+comments          +libcall           -python            +visual
+conceal           +linebreak         +python3           +visualextra
+cryptv            +lispindent        +quickfix          +viminfo
+cscope            +listcmds          +reltime           +vreplace
+cursorbind        +localmap          +rightleft         +wildignore
+cursorshape       +lua               +ruby              +wildmenu
+dialog_con        +menu              +scrollbind        +windows
+diff              +mksession         +signs             +writebackup
+digraphs          +modify_fname      +smartindent       -X11
-dnd               +mouse             -sound             -xfontset
-ebcdic            -mouseshape        +spell             -xim
+emacs_tags        +mouse_dec         +startuptime       -xpm
+eval              -mouse_gpm         +statusline        -xsmp
+ex_extra          -mouse_jsbterm     -sun_workshop      -xterm_clipboard
+extra_search      +mouse_netterm     +syntax            -xterm_save
   file vimrc di sistema: "$VIM/vimrc"
       file vimrc utente: "$HOME/.vimrc"
    II file vimrc utente: "~/.vim/vimrc"
        file exrc utente: "$HOME/.exrc"
        file dei default: "$VIMRUNTIME/defaults.vim"
         $VIM di riserva: "/usr/local/share/vim"
Compilazione: clang -c -I. -Iproto -DHAVE_CONFIG_H -DMACOS_X -DMACOS_X_DARWIN -g -O2 -U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=1
Link: clang -L. -fstack-protector-strong -L/usr/local/lib -L/usr/local/opt/libyaml/lib -L/usr/local/opt/openssl@1.1/lib -L/usr/local/opt/readline/lib -L/usr/local/lib -o vim -lm -lncurses -liconv -lintl -framework AppKit -L/usr/local/opt/lua/lib -llua5.4 -mmacosx-version-min=11.2 -fstack-protector-strong -L/usr/local/lib -L/usr/local/Cellar/perl/5.32.1_1/lib/perl5/5.32.1/darwin-thread-multi-2level/CORE -lperl -L/usr/local/opt/python@3.9/Frameworks/Python.framework/Versions/3.9/lib/python3.9/config-3.9-darwin -lpython3.9 -framework CoreFoundation -lruby.3.0 -L/usr/local/Cellar/ruby/3.0.1/lib

Output of YcmDebugInfo

Printing YouCompleteMe debug information...
-- Resolve completions: On demand
-- Client logfile: /var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/ycm_v053otmq.log
-- Server Python interpreter: /usr/local/opt/python@3.9/bin/python3.9
-- Server Python version: 3.9.5
-- Server has Clang support compiled in: False
-- Clang version: None
-- No extra configuration file found
-- GenericLSP completer debug information:
--   latexCompleter running
--   latexCompleter process ID: 89666
--   latexCompleter executable: ['/usr/local/bin/texlab']
--   latexCompleter logfiles:
--     /var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/latexcompleter_stderr17hqd1rv.log
--   latexCompleter Server State: Initialized
--   latexCompleter Project Directory: /Users/demichel/PAPERS/ISOBARICMD/LATEX
--   latexCompleter Settings: {}
-- Server running at: http://127.0.0.1:64159
-- Server process ID: 89665
-- Server logfiles:
--   /var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/ycmd_64159_stdout_eajav69t.log
--   /var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/ycmd_64159_stderr_vqmzpuy6.log

Output of YcmDiags

No warnings or errors detected.

Output of git rev-parse HEAD in YouCompleteMe installation directory

7c4d05375a09a871f618f9688c7af517d4e69b76

Contents of YCM, ycmd and completion engine logfiles

Reproduce the issue with vim -Nu /path/to/YCM/vimrc_ycm_minimal, which enabled debug logging and other useful diagnostics. Include a link to a [gist][] containing all of the log files listed by :YcmToggleLogs. file /var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/latexcompleter_stderr17hqd1rv.log:

---------------- file ```/var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/ycm_ufdu75sj.log```: ``` b'{"filepath": "/Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex", "line_num": 1, "column_num": 1, "working_dir": "/Users/demichel/PAPERS/ISOBARICMD/LATEX", "file_data": {"/Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex": {"contents": "\\\\documentclass[aps,pre,superscriptaddress]{revtex4-1} % chktex 8\\n%chktex-file 23\\n\\n\\\\usepackage{xpatch}\\n\\\\hbadness=10001\\n\\n\\\\usepackage{epsfig}\\n\\\\usepackage{graphicx}\\n\\\\usepackage{amssymb}\\n\\\\usepackage{amsmath}\\n\\\\usepackage{color}\\n\\\\usepackage{epstopdf}\\n\\\\usepackage{dsfont}\\n\\\\usepackage{romannum}\\n\\\\begin{document}\\n\\\\title{Isothermal-isobaric ensemble using the molecular virial: a reversible molecular dynamics implementation done right}\\n\\n\\\\author{Cristiano~De~Michele}\\n\\\\email{cristiano.demichele@uniroma1.it}\\n\\\\affiliation{Dipartimento di Fisica, ``Sapienza\'\' Universit\\\\`a di Roma, P.le A. Moro 2, 00185 Roma, Italy}\\n\\\\def\\\\ontop#1#2{\\\\genfrac{}{}{0pt}{}{#1}{#2}}\\n\\\\def\\\\bo#1{\\\\mathbf#1}\\n\\\\def\\\\mc#1{\\\\mathcal#1}\\n\\\\def\\\\pd#1{\\\\frac{\\\\partial}{\\\\partial#1}}\\n\\\\def\\\\dtd{\\\\frac{\\\\Delta{t}}{2}}\\n\\\\def\\\\dtq{\\\\frac{\\\\Delta{t}}{4}}\\n\\\\def\\\\wtL{\\\\widetilde{L}}\\n\\\\begin{abstract}\\n\\nImplementation of molecular isobaric ensemble based on MTK equations with holonomic constraints. \\n\\\\end{abstract}\\n\\n\\n\\\\date{\\\\today}\\n\\n\\n\\\\maketitle\\n\\n%\\\\section{Introduction}\\n% XXX\\n\\n\\\\section{Introduction}\\n\\nThe implementation of a reversible isobaric ensemble in molecular dynamics simulations has a longstanding history. Starting from \\nthe early approaches proposed by Andersen, Nos\\u00e8 and afterwards by Hoover, the most accurate method is the one in~\\\\cite{Martyna1994}.\\nThe method discussed in~\\\\cite{Martyna1994} can be straightforwardly applied to molecular systems were atoms are not subjected to holonomic constraints and it has been also extended system with holonomic constraints~\\\\cite{tuckermanBook,YuROLL2010} where the pressure is calculated by accounting explicitly for constraint forces (atomic pressure).\\n\\n\\\\section{Methods}\\n\\\\subsection{Models}\\n\\n\\\\subsection{Molecular MTK equation for isobaric ensemble}\\nMTK equations based on molecular pressure read~\\\\cite{tuckermanBook,Martyna1994,Martyna1996}:\\n\\\\begin{eqnarray}\\n \\\\bo{r}_{i,\\\\alpha} &=& \\\\frac{\\\\bo{p}_{i,\\\\alpha}}{m_{i,\\\\alpha}} + \\\\frac{p_\\\\epsilon}{W} \\\\bo{R}_i \\\\nonumber\\\\\\\\\\n \\\\dot{\\\\bo{p}}_{i,\\\\alpha} &=& \\\\bo{f}_{i,\\\\alpha} - \\\\left(1 + \\\\frac{3}{N_{fm}}\\\\right) \\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i} \\\\bo{P}_i -\\n \\\\frac{p_{\\\\eta_1}}{Q_{\\\\eta,1}} \\\\bo{p}_{i,\\\\alpha} \\\\nonumber\\\\\\\\\\n \\\\dot{V} &=& \\\\frac{3 V p_\\\\epsilon}{W}\\\\nonumber\\\\\\\\\\n \\\\dot{p}_\\\\epsilon &=& 3 V ({\\\\cal P}_{mol} - P ) + \\\\frac{3}{N_{fm}} \\\\sum_{i=1}^N \\\\frac{\\\\bo{P}^2_i}{M_i} - \\\\frac{p_{\\\\xi_1}}{Q_{\\\\xi,1}} p_\\\\epsilon\\\\nonumber \\\\\\\\\\n \\\\dot\\\\eta_j &=& \\\\frac{p_{\\\\eta_j}}{Q_{\\\\eta,j}} \\\\hspace{1.0cm} j=1,\\\\ldots,M_{\\\\eta}\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\eta_1} &=& \\\\left[ \\\\sum_{i,\\\\alpha} \\\\frac{\\\\bo{p}_i^2}{m_{i,\\\\alpha}} - N_f k_B T \\\\right] - \\\\frac{p_{\\\\eta_2}}{Q_{\\\\eta,2}} p_{\\\\eta_1}\\\\nonumber \\\\\\\\\\n \\\\dot{p}_{\\\\eta_j} &=& \\\\left[ \\\\frac{p^2_{\\\\eta_{j-1}}}{ Q_{\\\\eta,j-1}} - k_B T \\\\right ] - \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} p_{\\\\eta_j} \\\\hspace{0.5cm} j=2,\\\\ldots, M_\\\\eta-1\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\eta_M} &=& \\\\left [ \\\\frac{p^2_{\\\\eta_{M_\\\\eta-1}}}{Q_{\\\\eta,M_{\\\\eta}-1}} - k_B T \\\\right ]\\\\nonumber\\\\\\\\\\n \\\\dot\\\\xi_j &=& \\\\frac{p_{\\\\xi_j}}{Q_{\\\\xi,j}} \\\\hspace{1cm} j=1,\\\\ldots,M_{\\\\xi}\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\xi_1} &=& \\\\left[ \\\\frac{p_\\\\epsilon^2}{W} - k_B T \\\\right] - \\\\frac{p_{\\\\xi_2}}{Q_{\\\\xi,2}} p_{\\\\xi_1}\\\\nonumber \\\\\\\\\\n \\\\dot{p}_{\\\\xi_j} &=& \\\\left[ \\\\frac{p^2_{\\\\xi_{j-1}}}{Q_{\\\\xi,j-1}} - k_B T \\\\right ] - \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} p_{\\\\xi_j} \\\\hspace{0.5cm} j=2,\\\\ldots, M_\\\\xi-1\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\xi_{M_\\\\xi}} &=& \\\\left [ \\\\frac{p^2_{\\\\eta_{M_\\\\xi-1}}}{Q_{\\\\xi,M_{\\\\xi}-1}} - k_B T \\\\right ]\\n \\\\label{eq:MTK}\\n\\\\end{eqnarray}\\nwhere $i=1\\\\ldots N$ runs over all $N$ molecules and $\\\\alpha$ over all $n_i$ atoms within a molecule, $\\\\mathcal{P}_{mol}$ is the molecular pressure~\\\\cite{Ciccotti2004}, $\\\\bo{P}_i = \\\\sum_{\\\\alpha} \\\\bo{p}_{i,\\\\alpha}$,\\n\\\\begin{equation}\\n \\\\bo{R}_i = \\\\sum_{\\\\alpha} \\\\frac{m_{i,\\\\alpha} \\\\bo{r}_{i,\\\\alpha}}{M_i}\\n\\\\end{equation}\\nand we also coupled two independent Nos\\u00e8-Hoover chains to particles (whose associated variables are $\\\\eta_j$ and $p_{\\\\eta_j}$\\nwith $j=1\\\\ldots M_\\\\eta$) and to barostat (whose associated variables are $\\\\xi_j$ and $p_{\\\\xi_j}$ with $j=1\\\\ldots M_\\\\xi$).\\nIn Eq.~\\\\eqref{eq:MTK} $N_f$ is the number of degree of freedom in the system, $N_{fm}$ is the number of degree of freedom associated to molecules\\n(e.g.\\\\ if the total momentum is conserved $N_{fm}=3N - 3$), $Q_{\\\\eta,j}$ and $Q_{\\\\xi,j}$ are thermostat masses and\\n$W$ is the barostat mass.\\nTo numerically integrate these equations we implement a reversible integrator algorithm. If we define the state of the system as:\\n\\\\begin{equation}\\n \\\\Gamma = ( \\\\bo{r}^{N_{at}}, \\\\bo{p}^{N_{at}}, V, p_\\\\epsilon, \\\\eta_1\\\\ldots \\\\eta_{M_\\\\eta}, p_{\\\\eta_1}\\\\ldots p_{\\\\eta_{M_\\\\eta}},\\\\xi_1\\\\ldots \\\\xi_{M_\\\\xi}, p_{\\\\xi_1}\\\\ldots p_{\\\\xi_{M_\\\\xi}})\\n \\\\label{eq:Gamma}\\n\\\\end{equation}\\nwhere $N_{at}=\\\\sum_i n_i$ and $\\\\Gamma$ encompasses all the $6 N_{at} + 2 M_\\\\eta + 2 M_\\\\xi + 2 $ variables of the extended phase space,\\nthe equations of motion of the system can be written as follows:\\n\\\\begin{equation}\\n \\\\dot\\\\Gamma = i \\\\hat{L}\\\\, \\\\Gamma\\n\\\\label{eq:motion} \\n\\\\end{equation}\\nwhere $\\\\hat{L}$ is the Liouville operator, which can be explicitly written as follows:\\n\\\\begin{equation}\\n i\\\\hat{L} = i L_{NHB} + i L_{NHP} + i L_{p_\\\\epsilon} + i L_{{\\\\mathbf{p}},MTK} + i L_{\\\\mathbf p} + \\n i L_{\\\\mathbf{q},MTK} + i L_{\\\\mathbf{q}} + iL_{\\\\epsilon} \\n \\\\label{eq:liouville}\\n\\\\end{equation}\\nwhere $i L_{NHB}$ and $i L_{NHP}$ are the time evolution operators of the two Nos\\u00e8-Hoover chains~\\\\cite{tuckermanBook} coupled to \\nthe barostat and particles respectively, i.e\\n\\\\begin{eqnarray}\\n i L_{NHP} &=& - \\\\sum_{i=1}^{N_{at}} \\\\frac{p_{\\\\eta_1}}{Q_{\\\\eta,1}} \\\\bo{p}_i \\\\cdot \\\\pd{\\\\bo{p}_i} + \\\\sum_{j=1}^{M_\\\\eta} \\n \\\\frac{p_{\\\\eta_j}}{Q_{\\\\eta,j}} \\\\pd{p_{\\\\eta_j}} + \\\\sum_{j=1}^{M_{\\\\eta}-1} \\\\left ( G_{\\\\eta,j} - p_{\\\\eta_j}\\n \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\right ) \\\\pd{p_{\\\\eta_j}} + G_{\\\\eta,M} \\\\pd{p_{\\\\eta_{M_\\\\eta}}} \\\\label{eq:NHPop}\\\\\\\\ \\n i L_{NHB} &=& - \\\\frac{p_{\\\\xi_1}}{Q_{\\\\xi,1}} p_\\\\epsilon \\\\pd{p_\\\\epsilon} + \\\\sum_{j=1}^{M_\\\\xi} \\n \\\\frac{p_{\\\\xi_j}}{Q_{\\\\xi,j}} \\\\pd{p_{\\\\xi_j}} + \\\\sum_{j=1}^{M_{\\\\xi}-1} \\\\left ( G_{\\\\xi,j} - p_{\\\\xi_j}\\n \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\right ) \\\\pd{p_{\\\\xi_j}} + G_{\\\\xi,M} \\\\pd{p_{\\\\xi_{M_\\\\xi}}}\\n \\\\label{eq:NHBop}\\n\\\\end{eqnarray}\\nwhere \\n\\\\begin{eqnarray}\\n G_{\\\\eta,1} &=& \\\\sum_{i=1}^{N_{at}} \\\\frac{\\\\bo{p}_i^2}{m_i} - N_{f} k_B T\\\\\\\\\\n G_{\\\\eta,j}&=& \\\\frac{p_{\\\\eta_{j-1}}}{Q_{\\\\eta,j-1}} - k_B T \\\\\\\\\\n G_{\\\\xi,1} &=& \\\\frac{p_\\\\epsilon^2}{W} - k_B T\\\\\\\\\\n G_{\\\\xi,j}&=& \\\\frac{p_{\\\\xi_{j-1}}}{Q_{\\\\xi,j-1}} - k_B T \\n \\\\label{eq:Gdef}\\n\\\\end{eqnarray}\\nand:\\n\\\\begin{eqnarray}\\n i L_{\\\\bo{p},MTK} &=&\\\\sum_{i,\\\\alpha} \\\\left ( 1 + \\\\frac{3}{N_{fm}} \\\\right ) \\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i} \\\\bo{P}_i \\\\pd{\\\\bo{p}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\ \\n i L_{p_\\\\epsilon} &=& \\\\left [ 3 V ( \\\\mathcal{P}_{mol} - P ) + \\\\frac{3}{N_{fm}} \\\\sum_{i} \\\\frac{\\\\bo{P}_i^2}{M_i} \\\\right ] \\\\pd{p_\\\\epsilon}\\\\nonumber\\\\\\\\\\n i L_{\\\\bo{p}} &=& \\\\sum_{i,\\\\alpha} \\\\bo{f}_{i,\\\\alpha} \\\\pd{\\\\bo{p}_{i,\\\\alpha}}\\\\nonumber\\\\\\\\\\n i L_{\\\\bo{q},MTK} &=& \\\\sum_{i,\\\\alpha} \\\\frac{p_\\\\epsilon}{W} \\\\bo{R}_i \\\\pd{\\\\bo{r}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\\\n i L_{\\\\bo{q}} &=& \\\\frac{\\\\bo{p}_{i,\\\\alpha}}{m_{i,\\\\alpha}} \\\\pd{\\\\bo{r}_{i,\\\\alpha}}\\\\nonumber\\\\\\\\\\n i L_{\\\\epsilon} &=& \\\\frac{p_\\\\epsilon}{W} 3 V \\\\pd{V}\\n\\\\end{eqnarray}\\n\\\\section{Integration of the equations of motion}\\nTo integrate the equations of motion we apply the following decomposition of the operator $e^{i \\\\hat{L}}$:\\n\\\\begin{eqnarray}\\n e^{i \\\\hat{L}} &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{p},MTK} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t}\\n e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}}\\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}}\\\\nonumber\\\\\\\\ \\n && e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n \\\\label{eq:trotdec} \\n\\\\end{eqnarray}\\nwhere all the operators are applied from left to right. Latter decomposition is based on the general approach which we briefly illustrate\\nin the following. If $iL_g$ a generic operator which can be split as follows:\\n\\\\begin{eqnarray}\\n i L_g = i L_1 + i L_2 + i L_3\\n\\\\end{eqnarray}\\none can use a first Trotter decomposition of the time evolution operator $e^{i L_g \\\\Delta t}$ as follows:\\n\\n\\\\begin{equation}\\n e^{iL_g \\\\Delta t } \\\\approx e^{i L_3 \\\\frac{\\\\Delta t}{2}} e^{(iL_1 + iL_2) \\\\Delta t} e^{i L_3 \\\\frac{\\\\Delta t}{2} } \\n\\\\end{equation}\\n\\nand one can employ a further Trotter decomposition for the operator \\n $e^{(iL_1 + iL_2) \\\\Delta t}$, i.e.\\n \\\\begin{eqnarray}\\n e^{(iL_1 + iL_2) \\\\Delta t} \\\\approx e^{iL_2 \\\\frac{\\\\Delta t}{2}} e^{i L_1 \\\\Delta t} e^{i L_2 \\\\dtd}\\n \\\\end{eqnarray}\\n\\nso that we have for the operator $e^{i L_g \\\\Delta t}$:\\n\\n\\\\begin{equation}\\n e^{iL_g \\\\Delta t} \\\\approx e^{iL_3 \\\\dtd} e^{iL_2 d\\\\dtd} e^{i L_1 \\\\Delta t} e^{i L_2 \\\\dtd} e^{i L_3 \\\\dtd} \\n \\\\label{eq:trotterMF}\\n\\\\end{equation}\\nThe operators $e^{i L_{NHP} \\\\Delta t/2}$ and $ e^{i L_{NHB} \\\\Delta t/2} $ can be factorized by leveraging the Trotter theorem, i.e.~if \\nwe define the following two operators:\\n\\\\begin{eqnarray}\\n S_{NHP}(\\\\delta) &=& e^{ \\\\frac{\\\\delta}{2} G_{\\\\eta,M_\\\\eta} \\\\pd{p_{\\\\eta_{M_\\\\eta}}} } \\n \\\\prod_{j=M_{\\\\eta}}^{1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{\\\\eta_j} \\\\pd{p_{\\\\eta_j}}} e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n \\\\right ) \\\\prod_{i=1}^{N_{at}} e^{\\\\delta \\\\frac{p_{\\\\eta_1}}{Q_{\\\\eta,1}} \\\\bo{p}_{i}\\\\cdot \\\\pd{\\\\bo{p}_i} }\\\\nonumber\\\\\\\\\\n && \\\\prod_{j=1}^{M_{\\\\eta}} e^{- \\\\delta\\\\frac{p_{\\\\eta_j}}{Q_{\\\\eta,j}} \\\\pd{\\\\eta_j} }\\n \\\\prod^{M_{\\\\eta}}_{j=1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{_j}\\\\pd{p_{\\\\eta_j}}} \\n e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n \\\\right ) e^{ \\\\frac{\\\\delta}{2} G_{\\\\eta,M_\\\\eta} \\\\pd{p_{\\\\eta_{M_\\\\eta}}} }\\\\\\\\\\n S_{NHB}(\\\\delta) &=& e^{ \\\\frac{\\\\delta}{2} G_{,M_\\\\xi} \\\\pd{p_{_{M_\\\\xi}}} } \\n \\\\prod_{j=M_{\\\\xi}}^{1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{\\\\xi_j} \\\\pd{p_{\\\\xi_j}}} e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n \\\\right ) e^{\\\\delta \\\\frac{p_{\\\\xi_1}}{Q_{\\\\xi,1}} p_{\\\\epsilon} \\\\pd{p_\\\\epsilon} }\\\\nonumber\\\\\\\\\\n && \\\\prod_{j=1}^{M_{\\\\xi}} e^{-\\\\delta \\\\frac{p_{\\\\xi_j}}{Q_{\\\\xi,j}} \\\\pd{\\\\xi_j} }\\n \\\\prod^{M_{\\\\xi}}_{j=1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{\\\\xi_j}\\\\pd{p_{\\\\xi_j}}} \\n e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n \\\\right ) e^{ \\\\frac{\\\\delta}{2} G_{\\\\xi,{M_\\\\xi}} \\\\pd{p_{\\\\xi_{M_\\\\xi}}} }\\n \\\\label{eq:SNH}\\n\\\\end{eqnarray}\\none has:\\n\\\\begin{eqnarray}\\n e^{i L_{NHP} \\\\Delta t/2} &=& S_{NHP}(\\\\Delta t/2) + \\\\mathcal{O}({\\\\Delta t}^3) \\\\label{eq:trdecomNHP} \\\\\\\\ \\n e^{i L_{NHB} \\\\Delta t/2} &=& S_{NHB}(\\\\Delta t/2) + \\\\mathcal{O}({\\\\Delta t}^3)\\n \\\\label{eq:trdecomNHB}\\n\\\\end{eqnarray}\\n\\n$S_{NHP}$ and $S_{NHB}$ are called primitive decompositions of the operators $e^{i L_{NHP} \\\\Delta t/2}$ \\nand $ e^{i L_{NHB} \\\\Delta t/2} $ and they can be obtained by nothing that if one has a generic operator $i L_g$, \\nsuch that\\n\\\\begin{equation}\\n i \\\\wtL_g = i \\\\wtL_1 + i \\\\wtL_2 + i \\\\wtL_3 + i \\\\wtL_4\\n\\\\end{equation}\\nresorting to the same Trotter factorization used to arrive at Eq.~\\\\eqref{eq:trotterMF}\\nwith $iL_1 = i \\\\wtL$, $iL_2 = i \\\\wtL_2 + i\\\\wtL_3 $ and $ i L_3 = i \\\\wtL_4$, one obtains:\\n\\\\begin{eqnarray}\\n e^{i \\\\wtL_g \\\\Delta t} \\\\approx e^{i \\\\wtL_4 \\\\dtd} e^{ i (\\\\wtL_2 + i \\\\wtL_3) \\\\dtd} e^{ i \\\\wtL_1 \\\\Delta t}\\n e^{ i (\\\\wtL_2 + i \\\\wtL_3) \\\\dtd} e^{i \\\\wtL_4 \\\\dtd}\\n\\\\end{eqnarray}\\nand by factorizing again the operators $e^{ i (\\\\wtL_2 + i \\\\wtL_3) \\\\dtd}$, i.e.\\n\\\\begin{equation}\\n e^{(i\\\\wtL_2 + i \\\\wtL_3) \\\\dtd} \\\\approx e^{i \\\\wtL_3 \\\\dtq } e^{i \\\\wtL_2 \\\\dtd} e^{i \\\\wtL_3 \\\\dtq} \\n\\\\end{equation}\\none finally has:\\n\\\\begin{equation}\\n e^{i \\\\wtL_g \\\\Delta t} \\\\approx e^{i \\\\wtL_4 \\\\dtd} e^{i \\\\wtL_3 \\\\dtq } e^{i \\\\wtL_2 \\\\dtd} e^{i \\\\wtL_3 \\\\dtq} \\n e^{i \\\\wtL_4 \\\\dtd}\\n\\\\end{equation}\\n\\n% In Eqs.~\\\\eqref{eq:SYdecomNHP} and \\\\eqref{eq:SYdecomNHB} the weights $w_k$ can be chosen in such\\n% a way that an higher order integration scheme can be obtained. For example for a fourth-order scheme $n_{sy}=3$ and \\n% \\\\begin{eqnarray}\\n% w_1 &=& w_3 = \\\\frac{1}{2 - 2^{1/3}}\\\\nonumber\\\\\\\\\\n% w_2 &=& 1 - (w_1+w_3)\\n% \\\\label{eq:syweights}\\n% \\\\end{eqnarray}\\n% For a sixth-order scheme $n_{sy}=7$ and the weights are specified numerically~\\\\cite{tuckermanBook,Yoshida1990}. \\n% Note that if $n_{sy}=1$ and $w_1=1$ one recover the usual Trotter factorization.\\n\\nIn a similar fashion we can treat the operators $i L_{\\\\bo{q},MTK}$ and $i L_{\\\\bo{p},MTK}$. \\nIndeed, if we factorize the latter two operators as follows:\\n\\\\begin{eqnarray}\\n L_{\\\\bo{q},MTK} &=& \\\\sum_i L_{\\\\bo{q},MTK}^{(i)} = \\\\sum_i L_{\\\\bo{q},MTK0}^{(i)} + L_{\\\\bo{q},MTK1}^{(i)} \\\\nonumber\\\\\\\\ \\n L_{\\\\bo{p},MTK} &=& \\\\sum_i L_{\\\\bo{p},MTK}^{(i)} = \\\\sum_i L_{\\\\bo{p},MTK0}^{(i)} + L_{\\\\bo{p},MTK1}^{(i)} \\n \\\\label{eq:MTKop}\\n\\\\end{eqnarray}\\nwhere:\\n\\\\begin{eqnarray}\\n i L_{\\\\bo{q},MTK0}^{(i)} &=& \\\\sum_\\\\alpha \\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i}\\\\bo{r}_{i,\\\\alpha} \\\\pd{\\\\bo{r}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\ \\n i L_{\\\\bo{q},MTK1}^{(i)} &=& \\\\frac{p_\\\\epsilon}{W} \\\\sum_{\\\\ontop{\\\\alpha,\\\\alpha\'}{\\\\alpha\'\\\\neq\\\\alpha}}\\\\frac{m_{i,\\\\alpha\'}\\\\bo{r}_{i,\\\\alpha\'}}{M_i} \\\\pd{\\\\bo{r}_{i,\\\\alpha}} \\n \\\\label{eqMTKqsplit}\\n\\\\end{eqnarray}\\nand\\n\\\\begin{eqnarray}\\n i L_{\\\\bo{p},MTK0}^{(i)} &=& \\\\sum_{\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i}\\\\bo{p}_{i,\\\\alpha} \\\\pd{\\\\bo{p}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\ \\n i L_{\\\\bo{p},MTK1}^{(i)} &=& \\\\frac{p_\\\\epsilon}{W} \\\\sum_{\\\\ontop{\\\\alpha,\\\\alpha\'}{\\\\alpha\'\\\\neq\\\\alpha}}\\\\frac{m_{i,\\\\alpha}\\\\bo{p}_{i,\\\\alpha\'}}{M_i} \\\\pd{\\\\bo{p}_{i,\\\\alpha}} \\n \\\\label{eqMTKpsplit}\\n\\\\end{eqnarray}\\nagain we can employ the Trotter factorization to have:\\n\\\\begin{eqnarray}\\n e^{i L_{\\\\bo{q},MTK}^{(i)} \\\\Delta t/2} &=& S^{(i)}_{\\\\bo{q},MTK}(\\\\Delta t/2) + \\n \\\\mathcal{O}({\\\\Delta t}^3)\\\\label{eq:trdecomq} \\\\\\\\ \\n e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2} &=& S^{(i)}_{\\\\bo{p},MTK}(\\\\Delta t/2) + \\n \\\\mathcal{O}({\\\\Delta t}^3) \\n \\\\label{eq:trdecomp}\\n\\\\end{eqnarray}\\nwhere $S_{\\\\bo{q},MTK}(\\\\delta)$ and $S_{\\\\bo{p},MTK}(\\\\delta)$ are the following primitive factorizations of the operators \\n$e^{i L_{\\\\bo{q}}^{(i)}\\\\Delta t/2}$ and $e^{i L_{\\\\bo{p}}^{(i)} \\\\Delta t/2}$ respectively: \\n\\\\begin{eqnarray}\\n S^{(i)}_{\\\\bo{q},MTK}(\\\\delta) &=& \\\\prod_{\\\\alpha=1}^{n_i} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha\'}\\\\bo{r}_{i,\\\\alpha\'}}{M_i} \\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{r}_{i,\\\\alpha}}} \\n e^{i L_{\\\\bo{q},MTK0}^{(i)}\\\\delta } \\n \\\\prod_{\\\\alpha=n_i}^{1} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W}\\\\frac{m_{i,\\\\alpha\'}\\\\bo{r}_{i,\\\\alpha\'}}{M_i}\\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{r}_{i,\\\\alpha}}} \\n \\\\nonumber\\\\\\\\\\n S^{(i)}_{\\\\bo{p},MTK}(\\\\delta) &=& \\\\prod_{\\\\alpha=1}^{n_i} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}\\\\bo{p}_{i,\\\\alpha\'}}{M_i} \\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{p}_{i,\\\\alpha}}} \\n e^{i L_{\\\\bo{p},MTK0}^{(i)} \\\\delta} \\n \\\\prod_{\\\\alpha=n_i}^{1} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}\\\\bo{p}_{i,\\\\alpha\'}}{M_i} \\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{p}_{i,\\\\alpha}}}\\n \\\\label{eq:primfact} \\n\\\\end{eqnarray}\\n\\n\\\\noindent If the atoms of the $N$ molecules are subjected to the holonomic constraints:\\n\\\\begin{equation}\\n \\\\sigma_{i,\\\\beta}\\\\left (\\\\bo{r}_{i,1}\\\\ldots \\\\bo{r}_{i,n_i}\\\\right ) = 0\\n\\\\label{eq:constaints}\\n\\\\end{equation}\\nwhere $\\\\beta=1\\\\ldots l_i$ runs over all $l_i$ constraints of $i$-th molecule, the two operators $e^{i L_{\\\\bo{p}} \\\\Delta/2}$ \\nin Eq.~\\\\eqref{eq:trotdec} have to be modified to account for constraint forces\\nacting on each atom $(i,\\\\alpha)$ at time $t$ and $t+\\\\Delta t$. We assume that these forces can be expressed in terms of two different \\nsets of lagrangian multipliers ${\\\\lambda_{i,\\\\beta}(t)}$ and ${\\\\lambda_{i,\\\\beta}(t+\\\\Delta t)}$.\\nThe first operator from left in Eq.~\\\\eqref{eq:trotdec} evolves momenta according to forces evaluated at time $t$ and thus it becomes:\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{p}} \\\\Delta t/2} \\\\rightarrow e^{i L^{A}_{\\\\bo{p}} \\\\Delta t/2} = e^{\\\\left[\\\\sum_{i,\\\\alpha} \\\\bo{f}_{i,\\\\alpha} + \\n \\\\bo{f}^c_{i,\\\\alpha}(t) \\\\right] \\n \\\\pd{\\\\bo{p}_{i,\\\\alpha}}} \\n \\\\label{eq:op1constr}\\n\\\\end{equation}\\nwhere $\\\\bo{f}^c_{i,\\\\alpha}(t)$ is the sum of all constraint forces acting on the atom $(i,\\\\alpha)$, i.e. \\n\\\\begin{equation}\\n \\\\bo{f}^c_{i,\\\\alpha}(t) = \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t) \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{i,\\\\beta}(t)\\n\\\\end{equation}\\nand the second operator from left in Eq.~\\\\eqref{eq:trotdec} is modified as follows:\\n\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{p}} \\\\Delta t/2} \\\\rightarrow e^{i L^{B}_{\\\\bo{p}} \\\\Delta t/2} = e^{\\\\left[\\\\sum_{i,\\\\alpha} \\\\bo{f}_{i,\\\\alpha} \\n + \\\\bo{f}^c_{i,\\\\alpha}(t+\\\\Delta t) \\\\right]\\n \\\\pd{\\\\bo{p}_{i,\\\\alpha}}} \\n \\\\label{eq:op2constr}\\n\\\\end{equation}\\nwhere the constraint forces at time $t+\\\\Delta t$ are:\\n\\\\begin{equation}\\n \\\\bo{f}^c_{i,\\\\alpha}(t+\\\\Delta t) = \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t+\\\\Delta t) \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta (t+\\\\Delta t)\\n\\\\end{equation}\\nHence, the time evolution operator $e^{i \\\\hat{L}\\\\Delta t}$ becomes:\\n \\\\begin{eqnarray}\\n e^{i \\\\hat{L}} &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L^A_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t}\\n e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L^B_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}}\\\\nonumber\\\\\\\\ \\n && e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n \\\\label{eq:trotdec_constr}\\n\\\\end{eqnarray}\\n%We note that after the action of the operator $e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}}$ \\n%the correction of the position of an atom $(i,\\\\alpha)$ due to constraint forces \\n%depends only on constraint forces acting on $(i,\\\\alpha)$, i.e.\\\\@:\\n%\\\\begin{equation}\\n% e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\\\bo{r}_{i,\\\\alpha} = \\\\bo{r}_{i,\\\\alpha}^{nc} + \\\\frac{{\\\\Delta t}^2}{2} \\n% \\\\sum_\\\\beta \\\\tilde\\\\lambda_{i,\\\\beta}(t)\\n% \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{i,\\\\beta}(t)\\n%\\\\end{equation}\\n%where $\\\\bo{r}_{i,\\\\alpha}^{nc}$ is the position of the atom in absence of constraint forces and \\n%$\\\\tilde\\\\lambda_{i,\\\\beta}(t)$ are a set of unknown coefficients (see below on how to estimate them).\\n%Similarly, after the action of the operator $e^{i L_{\\\\bo{p},MTK} \\\\frac{\\\\Delta t}{2}}$ \\n%the correction of the velocity of an atom $(i,\\\\alpha)$ due to constraint forces \\n%depends only on constraint forces acting on $(i,\\\\alpha)$, i.e.\\\\@:\\n%\\\\begin{equation}\\n% e^{i L_{\\\\bo{p},MTK} \\\\frac{\\\\Delta t}{2}} \\\\bo{v}_{i,\\\\alpha} = \\\\bo{v}_{i,\\\\alpha}^{nc} + \\\\frac{\\\\Delta t}{2} \\\\sum_\\\\beta \\\\tilde\\\\lambda_{i,\\\\beta}(t+\\\\Delta t)\\n% \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{i,\\\\beta}(t)\\n%\\\\end{equation}\\n%where $\\\\bo{v}_{i,\\\\alpha}^{nc}$ is the velocity of the atom in absence of constraint forces.\\n\\nA convenient way to estimate the constraint forces is provided by \\nSHAKE algorithm~\\\\cite{Ciccotti1977,Ciccotti1982}. \\nWithin this approach, the full propagator is rewritten as follows:\\n\\\\begin{eqnarray}\\n e^{i \\\\hat{L}} &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t} \\\\hat{\\\\mathcal{S}}_{r}\\n e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\\\nonumber\\\\\\\\\\n && \\\\hat{\\\\mathcal{S}}_{v}\\n e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}} e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n \\\\label{eq:trotdecshk}\\n\\\\end{eqnarray}\\nwhere the operator $\\\\hat{\\\\mathcal{S}}_{r}$ adjust positions and velocities of atoms to take into \\naccount the action of the constraint forces at time $t$ and the operator $\\\\hat{\\\\mathcal{S}}_{v}$ adjust the velocities \\naccording to constraint forces at time $t+\\\\Delta t$. Their action is described in detail in the following.\\n\\nThe operator $\\\\hat{\\\\mathcal{S}}_{r}$ adjusts first the positions of particles to fulfill the constraints in Eq.~\\\\eqref{eq:constaints}, then\\nupdates the velocities accordingly, i.e.~if the displacement of particle $(i,\\\\alpha)$ is $\\\\delta \\\\bo{r}_{i,\\\\alpha}$ the corresponding \\nvariation of its velocity is $\\\\delta \\\\bo{r}_{i,\\\\alpha}/\\\\Delta t$.\\nThe adjustment of the positions is achieved iteratively as in the standard SHAKE approach~\\\\cite{Ciccotti1977}, where\\nthe bonds, which always involve pair of atoms, are adjusted one by one until convergence is reached. \\n%We note that the displacement of two atoms involved in a bond\\n%has to occur along a direction calculated according to the positions of particles prior to the application of the operator\\n%$e^{i L_{\\\\bo{q}} \\\\Delta t}$. \\n%The calculation of the constraint forces $\\\\bo{f}^c_{i,\\\\alpha}(t)$ for updating the velocities \\n%requires some care.\\nIf we define the states $\\\\Gamma $, $\\\\Gamma\'$ and $\\\\Gamma\'\'$ as follows:\\n\\\\begin{eqnarray}\\n \\\\Gamma\' &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}} e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\\\Gamma(t) \\\\\\\\\\n \\\\Gamma\'\' &=& e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t} \\\\, \\\\Gamma\'\\\\nonumber\\\\\\\\\\n %\\\\Gamma\'\' &=& e^{i L_{\\\\bo{q},MTK} \\\\Delta t /2}\\\\, \\\\Gamma\'\\n \\\\label{eq:rnoc}\\n\\\\end{eqnarray}\\nafter the application of $\\\\hat{\\\\mathcal{S}}_r$ the positions in state $\\\\Gamma\'\'$ will be adjusted\\nto fulfill all the holonomic constraints in Eq.~\\\\eqref{eq:constaints}, i.e.~if\\n\\\\begin{equation}\\n \\\\Gamma\'\'\' = \\\\hat{\\\\mathcal{S}}_{r} \\\\,\\\\Gamma\'\'% chktex 23 \\n \\\\label{eq:rp} \\n\\\\end{equation}\\nthe positions in state $\\\\Gamma\'\'\'$ fulfill the constraints in Eq.~\\\\eqref{eq:constaints}. % chktex 23\\n\\nLet\'s indicate with $\\\\bo{r}_{i,\\\\alpha}\'$, $\\\\bo{r}_{i,\\\\alpha}\'\'$ and $\\\\bo{r}_{i,\\\\alpha}\'\'\'$ \\nthe positions corresponding to states $\\\\Gamma\'$, $\\\\Gamma\'\'$ and $\\\\Gamma\'\'\'$ respectively.\\nThe operator $ \\\\hat{\\\\mathcal{S}}_{r} $ adjusts the positions by finding a set of lagrangian multipliers\\nsuch that:\\n\\\\begin{equation}\\n \\\\hat{\\\\mathcal{S}}_{r} \\\\bo{r}_{i,\\\\alpha}\'\' = \\\\bo{r}_{i,\\\\alpha}\'\' + \\\\frac{{\\\\Delta t}^2}{2} \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t) \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta (t)\\n \\\\label{eq:Sract}\\n\\\\end{equation}\\nwhere we note that $\\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta(t)$ must be calculated by using positions in state $\\\\Gamma\'$.\\nWe note that the operator $\\\\hat{S}_r$ can be safely applied before the operator $ e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} $,\\nsince latter one does not alter the bonds, indeed:\\n\\\\begin{equation}\\n %\\\\dot{\\\\bo{r}}_{i,\\\\alpha} - \\\\dot{\\\\bo{r}}_{i,\\\\alpha\'} = \\n i L_{\\\\bo{q}}^{(i),MTK} (\\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'}) = 0\\n\\\\end{equation}\\nhence\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{q},MTK}^{(i)} \\\\frac{\\\\Delta t}{2}} (\\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'} ) = \\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'}\\n \\\\label{eq:distinv}\\n\\\\end{equation}\\nand from Eq.~\\\\eqref{eq:SYdecomq} one has that:\\n\\\\begin{equation}\\n \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{q},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n (\\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'} ) = \\n \\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'} + \\\\mathcal{O}({\\\\Delta t}^3)\\n \\\\label{eq:rconstrconserv}\\n\\\\end{equation}\\n%We note that the quantities $\\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta(t)$ are calculated by using positions \\n%in state $\\\\Gamma_0$.\\n\\nAs in the case of operator $\\\\hat{S}_r$ the adjustments of the velocities through the operator $\\\\hat{\\\\mathcal{S}}_{v}$\\ncan be viewed as the result of the action of constraint forces $\\\\bo{f}^c_{i,\\\\alpha}(t+\\\\Delta t)$ as the ones in\\nthe operator $e^{i L^B_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}}$ in Eq.~\\\\eqref{eq:trotdec_constr}.\\nIf we define the states $\\\\Gamma^{\\\\Romannum{4}}$ and $\\\\Gamma^{\\\\Romannum{5}}$ as follows:\\n\\\\begin{eqnarray}\\n \\\\Gamma^{\\\\Romannum{4}} &=& e^{i L_{\\\\bo{q}, MTK} \\\\frac{\\\\Delta t}{2}} \\\\Gamma\'\'\'\\\\\\\\\\n \\\\Gamma^{\\\\Romannum{5}} &=&e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\\\Gamma^{\\\\Romannum{4}} \\n \\\\label{eq:gamma45}\\n\\\\end{eqnarray}\\nthe action of the operator $\\\\hat{S}_v$ can be expressed as:\\n\\\\begin{equation}\\n \\\\Gamma^{\\\\Romannum{6}} = \\\\hat{S}_v \\\\Gamma^{\\\\Romannum{5}} = \\\\bo{v}^{\\\\Romannum{5}}_{i,\\\\alpha} + \\\\frac{{\\\\Delta t}^2}{2} \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t+\\\\Delta t)\\n \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{\\\\beta}(t+\\\\Delta t)\\n\\\\end{equation}\\nwhere we note that $\\\\sigma_{\\\\beta}(t+\\\\Delta t)$ is calculated by using positions in state $\\\\Gamma^{\\\\Romannum{4}}$.\\nWe note that the operator $\\\\hat{S}_v$ can be safely applied after the operator $e^{i L_{\\\\bo{p}}^{(i)} \\\\Delta t/2}$, since latter\\none does not alter the relative velocity of atoms belonging to the same molecule, i.e.: \\n\\\\begin{equation}\\n i L_{\\\\bo{p},MTK}^{(i)} (\\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'}) = 0\\n\\\\end{equation}\\nso that\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\frac{\\\\Delta t}{2}} (\\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'} ) = \\n \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'}\\n \\\\label{eq:velinv}\\n\\\\end{equation}\\nand according to Eq.~\\\\eqref{eq:SYdecomp}:\\n\\\\begin{equation}\\n \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{p},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n (\\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\n \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'} ) = \\n \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'} + \\\\mathcal{O}({\\\\Delta t}^3)\\n \\\\label{eq:vconstrconserv}\\n\\\\end{equation}\\n\\n\\\\subsection{Suzuki-Yoshida decompositions}\\n\\nThe factorization of operators $e^{i L_{NHP}\\\\Delta t/2}$ and $e^{i L_{NHB}\\\\Delta t/2}$ in Eqs.~\\\\eqref{eq:trdecomNHP} and~\\\\eqref{eq:trdecomNHB}\\nand can be further improved by using a Suzuki-Yoshida (SY) decomposition and introducing RESPA \\nas discussed in~\\\\cite{tuckermanBook,Suzuki1991,Yoshida1990}.\\nSince $S_{NHP}$ and $S_{NHB}$ are primitive factorizations of the operators $e^{i L_{NHP}}$ and $e^{i L_{NHB}}$ we can \\nwrite~\\\\cite{tuckermanBook}:\\n%applying the SY decomposition to the operators $e^{i L_{NHP}\\\\Delta t/2}$ and $e^{i L_{NHB} \\\\Delta t/2}$ and by also introducing \\n%a RESPA with $n$ steps~\\\\cite{tuckermanBook} one has:\\n\\\\begin{eqnarray}\\ne^{i L_{NHP} \\\\Delta t/2} &\\\\approx& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S_{NHP}(w_k \\\\Delta t/2n)\\\\right]}^{n}\\\\label{eq:SYdecomNHP} \\\\\\\\ \\n e^{i L_{NHB} \\\\Delta t/2} &\\\\approx& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S_{NHB}(w_k \\\\Delta t/2n)\\\\right]}^{n} \\n \\\\label{eq:SYdecomNHB}\\n\\\\end{eqnarray}\\n\\nSimilarly, we can apply the Suzuki-Yoshida decomposition to the operators $e^{i L_{\\\\bo{q},MTK}^{(i)}\\\\Delta t/2}$ and \\n$e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2}$ and we can also introduce a RESPA with $n$ steps.\\nIn this case the primitive decompositions of the operators $e^{i L_{\\\\bo{q},MTK}^{(i)}\\\\Delta t/2}$ and \\n$e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2}$ are $S^{(i)}_{\\\\bo{q},MTK}$ and $S^{(i)}_{\\\\bo{p},MTK}$ so that:\\n\\n\\\\begin{eqnarray}\\n e^{i L_{\\\\bo{q},MTK}^{(i)} \\\\Delta t/2} &=& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{q},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n + \\n \\\\mathcal{O}({\\\\Delta t}^p) \\\\label{eq:SYdecomp} \\\\\\\\ \\n e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2} &\\\\approx& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{p},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n \\n + \\\\mathcal{O}({\\\\Delta t}^p)\\n \\\\label{eq:SYdecomq}\\n\\\\end{eqnarray}\\nwhere $p$ is the order of the SY decomposition.\\nFinally, we note that thanks to SY decomposition the error in Eqs.~\\\\eqref{eq:rconstrconserv} and~\\\\eqref{eq:vconstrconserv} \\nbecomes $\\\\mathcal{O}({\\\\Delta t}^p)$, i.e.~it significantly reduces.\\n\\n\\\\section{Results and Discussion}\\nXXX\\n\\n\\\\section{Conclusions}\\n\\nYYY\\n\\n%\\\\subsection{Acknowledgments}\\n%bla, bla\\\\ldots\\n\\\\bibliography{isobaric}\\n\\\\end{document}\\n", "filetypes": ["tex" ``` --------------------------- file ```/var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/ycmd_64159_stderr_vqmzpuy6.log```: ``` 2021-05-13 11:00:18,049 - DEBUG - No global extra conf, not calling method YcmCorePreload 2021-05-13 11:00:18,074 - INFO - Completion config: 50, detailing 10 candiates 2021-05-13 11:00:18,074 - INFO - Completion config: 50, detailing 10 candiates 2021-05-13 11:00:18,074 - INFO - Completion config: 50, detailing 10 candiates 2021-05-13 11:00:18,075 - INFO - Completion config: 50, detailing 10 candiates 127.0.0.1 - - [13/May/2021 11:00:18] "GET /ready HTTP/1.1" 200 4 2021-05-13 11:00:18,140 - INFO - Completion config: 50, detailing 10 candiates 127.0.0.1 - - [13/May/2021 11:00:18] "GET /signature_help_available?subserver=tex HTTP/1.1" 200 23 2021-05-13 11:00:18,142 - DEBUG - Event name: BufferVisit 127.0.0.1 - - [13/May/2021 11:00:18] "POST /event_notification HTTP/1.1" 200 2 2021-05-13 11:00:18,143 - DEBUG - Event name: FileReadyToParse 2021-05-13 11:00:18,143 - INFO - Adding buffer identifiers for file: /Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex 2021-05-13 11:00:18,149 - INFO - Starting latexCompleter: ['/usr/local/bin/texlab'] 2021-05-13 11:00:18,153 - INFO - latexCompleter started with PID 89666 2021-05-13 11:00:18,153 - DEBUG - TX: Sending message: b'Content-Length: 1261\r\n\r\n{"id":1,"jsonrpc":"2.0","method":"initialize","params":{"capabilities":{"textDocument":{"codeAction":{"codeActionLiteralSupport":{"codeActionKind":{"valueSet":["","quickfix","refactor","refactor.extract","refactor.inline","refactor.rewrite","source","source.organizeImports"]}}},"completion":{"completionItem":{"documentationFormat":["plaintext","markdown"]},"completionItemKind":{"valueSet":[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]}},"documentSymbol":{"hierarchicalDocumentSymbolSupport":false,"labelSupport":false,"symbolKind":{"valueSet":[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]}},"hover":{"contentFormat":["plaintext","markdown"]},"signatureHelp":{"signatureInformation":{"documentationFormat":["plaintext","markdown"],"parameterInformation":{"labelOffsetSupport":true}}},"synchronization":{"didSave":true}},"workspace":{"applyEdit":true,"didChangeWatchedFiles":{"dynamicRegistration":true},"symbol":{"symbolKind":{"valueSet":[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]}},"workspaceEdit":{"documentChanges":true}}},"initializationOptions":{},"processId":89665,"rootPath":"/Users/demichel/PAPERS/ISOBARICMD/LATEX","rootUri":"file:///Users/demichel/PAPERS/ISOBARICMD/LATEX"}}' 127.0.0.1 - - [13/May/2021 11:00:18] "POST /event_notification HTTP/1.1" 200 2 2021-05-13 11:00:18,234 - DEBUG - RX: Received message: b'{"jsonrpc":"2.0","result":{"capabilities":{"completionProvider":{"resolveProvider":true,"triggerCharacters":["\\\\","{","}","@","/"," "]},"definitionProvider":true,"documentFormattingProvider":true,"documentHighlightProvider":true,"documentLinkProvider":{"resolveProvider":false},"documentSymbolProvider":true,"foldingRangeProvider":true,"hoverProvider":true,"referencesProvider":true,"renameProvider":{"prepareProvider":true},"textDocumentSync":{"change":1,"openClose":true,"save":{"includeText":false}},"workspaceSymbolProvider":true},"serverInfo":{"name":"TexLab","version":"2.2.2"}},"id":1}' 2021-05-13 11:00:18,234 - INFO - None: Language server requires resolve request 2021-05-13 11:00:18,234 - INFO - None: Language server requires sync type of Full 2021-05-13 11:00:18,234 - DEBUG - latex: Server declares trigger characters: ['\\', '{', '}', '@', '/', ' '] 2021-05-13 11:00:18,234 - INFO - latex: Using trigger characters for semantic triggers: \,{,},@,/, 2021-05-13 11:00:18,234 - DEBUG - latex: Server declares signature trigger characters: [] 2021-05-13 11:00:18,235 - DEBUG - TX: Sending notification: b'Content-Length: 52\r\n\r\n{"jsonrpc":"2.0","method":"initialized","params":{}}' 2021-05-13 11:00:18,235 - DEBUG - TX: Sending notification: b'Content-Length: 86\r\n\r\n{"jsonrpc":"2.0","method":"workspace/didChangeConfiguration","params":{"settings":{}}}' 2021-05-13 11:00:18,235 - DEBUG - Refreshing file /Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex: State is Open/action Open 2021-05-13 11:00:18,235 - DEBUG - TX: Sending notification: b'Content-Length: 30168\r\n\r\n{"jsonrpc":"2.0","method":"textDocument/didOpen","params":{"textDocument":{"languageId":"tex","text":"\\\\documentclass[aps,pre,superscriptaddress]{revtex4-1} % chktex 8\\n%chktex-file 23\\n\\n\\\\usepackage{xpatch}\\n\\\\hbadness=10001\\n\\n\\\\usepackage{epsfig}\\n\\\\usepackage{graphicx}\\n\\\\usepackage{amssymb}\\n\\\\usepackage{amsmath}\\n\\\\usepackage{color}\\n\\\\usepackage{epstopdf}\\n\\\\usepackage{dsfont}\\n\\\\usepackage{romannum}\\n\\\\begin{document}\\n\\\\title{Isothermal-isobaric ensemble using the molecular virial: a reversible molecular dynamics implementation done right}\\n\\n\\\\author{Cristiano~De~Michele}\\n\\\\email{cristiano.demichele@uniroma1.it}\\n\\\\affiliation{Dipartimento di Fisica, ``Sapienza\'\' Universit\\\\`a di Roma, P.le A. Moro 2, 00185 Roma, Italy}\\n\\\\def\\\\ontop#1#2{\\\\genfrac{}{}{0pt}{}{#1}{#2}}\\n\\\\def\\\\bo#1{\\\\mathbf#1}\\n\\\\def\\\\mc#1{\\\\mathcal#1}\\n\\\\def\\\\pd#1{\\\\frac{\\\\partial}{\\\\partial#1}}\\n\\\\def\\\\dtd{\\\\frac{\\\\Delta{t}}{2}}\\n\\\\def\\\\dtq{\\\\frac{\\\\Delta{t}}{4}}\\n\\\\def\\\\wtL{\\\\widetilde{L}}\\n\\\\begin{abstract}\\n\\nImplementation of molecular isobaric ensemble based on MTK equations with holonomic constraints. \\n\\\\end{abstract}\\n\\n\\n\\\\date{\\\\today}\\n\\n\\n\\\\maketitle\\n\\n%\\\\section{Introduction}\\n% XXX\\n\\n\\\\section{Introduction}\\n\\nThe implementation of a reversible isobaric ensemble in molecular dynamics simulations has a longstanding history. Starting from \\nthe early approaches proposed by Andersen, Nos\\u00e8 and afterwards by Hoover, the most accurate method is the one in~\\\\cite{Martyna1994}.\\nThe method discussed in~\\\\cite{Martyna1994} can be straightforwardly applied to molecular systems were atoms are not subjected to holonomic constraints and it has been also extended system with holonomic constraints~\\\\cite{tuckermanBook,YuROLL2010} where the pressure is calculated by accounting explicitly for constraint forces (atomic pressure).\\n\\n\\\\section{Methods}\\n\\\\subsection{Models}\\n\\n\\\\subsection{Molecular MTK equation for isobaric ensemble}\\nMTK equations based on molecular pressure read~\\\\cite{tuckermanBook,Martyna1994,Martyna1996}:\\n\\\\begin{eqnarray}\\n \\\\bo{r}_{i,\\\\alpha} &=& \\\\frac{\\\\bo{p}_{i,\\\\alpha}}{m_{i,\\\\alpha}} + \\\\frac{p_\\\\epsilon}{W} \\\\bo{R}_i \\\\nonumber\\\\\\\\\\n \\\\dot{\\\\bo{p}}_{i,\\\\alpha} &=& \\\\bo{f}_{i,\\\\alpha} - \\\\left(1 + \\\\frac{3}{N_{fm}}\\\\right) \\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i} \\\\bo{P}_i -\\n \\\\frac{p_{\\\\eta_1}}{Q_{\\\\eta,1}} \\\\bo{p}_{i,\\\\alpha} \\\\nonumber\\\\\\\\\\n \\\\dot{V} &=& \\\\frac{3 V p_\\\\epsilon}{W}\\\\nonumber\\\\\\\\\\n \\\\dot{p}_\\\\epsilon &=& 3 V ({\\\\cal P}_{mol} - P ) + \\\\frac{3}{N_{fm}} \\\\sum_{i=1}^N \\\\frac{\\\\bo{P}^2_i}{M_i} - \\\\frac{p_{\\\\xi_1}}{Q_{\\\\xi,1}} p_\\\\epsilon\\\\nonumber \\\\\\\\\\n \\\\dot\\\\eta_j &=& \\\\frac{p_{\\\\eta_j}}{Q_{\\\\eta,j}} \\\\hspace{1.0cm} j=1,\\\\ldots,M_{\\\\eta}\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\eta_1} &=& \\\\left[ \\\\sum_{i,\\\\alpha} \\\\frac{\\\\bo{p}_i^2}{m_{i,\\\\alpha}} - N_f k_B T \\\\right] - \\\\frac{p_{\\\\eta_2}}{Q_{\\\\eta,2}} p_{\\\\eta_1}\\\\nonumber \\\\\\\\\\n \\\\dot{p}_{\\\\eta_j} &=& \\\\left[ \\\\frac{p^2_{\\\\eta_{j-1}}}{ Q_{\\\\eta,j-1}} - k_B T \\\\right ] - \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} p_{\\\\eta_j} \\\\hspace{0.5cm} j=2,\\\\ldots, M_\\\\eta-1\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\eta_M} &=& \\\\left [ \\\\frac{p^2_{\\\\eta_{M_\\\\eta-1}}}{Q_{\\\\eta,M_{\\\\eta}-1}} - k_B T \\\\right ]\\\\nonumber\\\\\\\\\\n \\\\dot\\\\xi_j &=& \\\\frac{p_{\\\\xi_j}}{Q_{\\\\xi,j}} \\\\hspace{1cm} j=1,\\\\ldots,M_{\\\\xi}\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\xi_1} &=& \\\\left[ \\\\frac{p_\\\\epsilon^2}{W} - k_B T \\\\right] - \\\\frac{p_{\\\\xi_2}}{Q_{\\\\xi,2}} p_{\\\\xi_1}\\\\nonumber \\\\\\\\\\n \\\\dot{p}_{\\\\xi_j} &=& \\\\left[ \\\\frac{p^2_{\\\\xi_{j-1}}}{Q_{\\\\xi,j-1}} - k_B T \\\\right ] - \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} p_{\\\\xi_j} \\\\hspace{0.5cm} j=2,\\\\ldots, M_\\\\xi-1\\\\nonumber\\\\\\\\\\n \\\\dot{p}_{\\\\xi_{M_\\\\xi}} &=& \\\\left [ \\\\frac{p^2_{\\\\eta_{M_\\\\xi-1}}}{Q_{\\\\xi,M_{\\\\xi}-1}} - k_B T \\\\right ]\\n \\\\label{eq:MTK}\\n\\\\end{eqnarray}\\nwhere $i=1\\\\ldots N$ runs over all $N$ molecules and $\\\\alpha$ over all $n_i$ atoms within a molecule, $\\\\mathcal{P}_{mol}$ is the molecular pressure~\\\\cite{Ciccotti2004}, $\\\\bo{P}_i = \\\\sum_{\\\\alpha} \\\\bo{p}_{i,\\\\alpha}$,\\n\\\\begin{equation}\\n \\\\bo{R}_i = \\\\sum_{\\\\alpha} \\\\frac{m_{i,\\\\alpha} \\\\bo{r}_{i,\\\\alpha}}{M_i}\\n\\\\end{equation}\\nand we also coupled two independent Nos\\u00e8-Hoover chains to particles (whose associated variables are $\\\\eta_j$ and $p_{\\\\eta_j}$\\nwith $j=1\\\\ldots M_\\\\eta$) and to barostat (whose associated variables are $\\\\xi_j$ and $p_{\\\\xi_j}$ with $j=1\\\\ldots M_\\\\xi$).\\nIn Eq.~\\\\eqref{eq:MTK} $N_f$ is the number of degree of freedom in the system, $N_{fm}$ is the number of degree of freedom associated to molecules\\n(e.g.\\\\ if the total momentum is conserved $N_{fm}=3N - 3$), $Q_{\\\\eta,j}$ and $Q_{\\\\xi,j}$ are thermostat masses and\\n$W$ is the barostat mass.\\nTo numerically integrate these equations we implement a reversible integrator algorithm. If we define the state of the system as:\\n\\\\begin{equation}\\n \\\\Gamma = ( \\\\bo{r}^{N_{at}}, \\\\bo{p}^{N_{at}}, V, p_\\\\epsilon, \\\\eta_1\\\\ldots \\\\eta_{M_\\\\eta}, p_{\\\\eta_1}\\\\ldots p_{\\\\eta_{M_\\\\eta}},\\\\xi_1\\\\ldots \\\\xi_{M_\\\\xi}, p_{\\\\xi_1}\\\\ldots p_{\\\\xi_{M_\\\\xi}})\\n \\\\label{eq:Gamma}\\n\\\\end{equation}\\nwhere $N_{at}=\\\\sum_i n_i$ and $\\\\Gamma$ encompasses all the $6 N_{at} + 2 M_\\\\eta + 2 M_\\\\xi + 2 $ variables of the extended phase space,\\nthe equations of motion of the system can be written as follows:\\n\\\\begin{equation}\\n \\\\dot\\\\Gamma = i \\\\hat{L}\\\\, \\\\Gamma\\n\\\\label{eq:motion} \\n\\\\end{equation}\\nwhere $\\\\hat{L}$ is the Liouville operator, which can be explicitly written as follows:\\n\\\\begin{equation}\\n i\\\\hat{L} = i L_{NHB} + i L_{NHP} + i L_{p_\\\\epsilon} + i L_{{\\\\mathbf{p}},MTK} + i L_{\\\\mathbf p} + \\n i L_{\\\\mathbf{q},MTK} + i L_{\\\\mathbf{q}} + iL_{\\\\epsilon} \\n \\\\label{eq:liouville}\\n\\\\end{equation}\\nwhere $i L_{NHB}$ and $i L_{NHP}$ are the time evolution operators of the two Nos\\u00e8-Hoover chains~\\\\cite{tuckermanBook} coupled to \\nthe barostat and particles respectively, i.e\\n\\\\begin{eqnarray}\\n i L_{NHP} &=& - \\\\sum_{i=1}^{N_{at}} \\\\frac{p_{\\\\eta_1}}{Q_{\\\\eta,1}} \\\\bo{p}_i \\\\cdot \\\\pd{\\\\bo{p}_i} + \\\\sum_{j=1}^{M_\\\\eta} \\n \\\\frac{p_{\\\\eta_j}}{Q_{\\\\eta,j}} \\\\pd{p_{\\\\eta_j}} + \\\\sum_{j=1}^{M_{\\\\eta}-1} \\\\left ( G_{\\\\eta,j} - p_{\\\\eta_j}\\n \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\right ) \\\\pd{p_{\\\\eta_j}} + G_{\\\\eta,M} \\\\pd{p_{\\\\eta_{M_\\\\eta}}} \\\\label{eq:NHPop}\\\\\\\\ \\n i L_{NHB} &=& - \\\\frac{p_{\\\\xi_1}}{Q_{\\\\xi,1}} p_\\\\epsilon \\\\pd{p_\\\\epsilon} + \\\\sum_{j=1}^{M_\\\\xi} \\n \\\\frac{p_{\\\\xi_j}}{Q_{\\\\xi,j}} \\\\pd{p_{\\\\xi_j}} + \\\\sum_{j=1}^{M_{\\\\xi}-1} \\\\left ( G_{\\\\xi,j} - p_{\\\\xi_j}\\n \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\right ) \\\\pd{p_{\\\\xi_j}} + G_{\\\\xi,M} \\\\pd{p_{\\\\xi_{M_\\\\xi}}}\\n \\\\label{eq:NHBop}\\n\\\\end{eqnarray}\\nwhere \\n\\\\begin{eqnarray}\\n G_{\\\\eta,1} &=& \\\\sum_{i=1}^{N_{at}} \\\\frac{\\\\bo{p}_i^2}{m_i} - N_{f} k_B T\\\\\\\\\\n G_{\\\\eta,j}&=& \\\\frac{p_{\\\\eta_{j-1}}}{Q_{\\\\eta,j-1}} - k_B T \\\\\\\\\\n G_{\\\\xi,1} &=& \\\\frac{p_\\\\epsilon^2}{W} - k_B T\\\\\\\\\\n G_{\\\\xi,j}&=& \\\\frac{p_{\\\\xi_{j-1}}}{Q_{\\\\xi,j-1}} - k_B T \\n \\\\label{eq:Gdef}\\n\\\\end{eqnarray}\\nand:\\n\\\\begin{eqnarray}\\n i L_{\\\\bo{p},MTK} &=&\\\\sum_{i,\\\\alpha} \\\\left ( 1 + \\\\frac{3}{N_{fm}} \\\\right ) \\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i} \\\\bo{P}_i \\\\pd{\\\\bo{p}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\ \\n i L_{p_\\\\epsilon} &=& \\\\left [ 3 V ( \\\\mathcal{P}_{mol} - P ) + \\\\frac{3}{N_{fm}} \\\\sum_{i} \\\\frac{\\\\bo{P}_i^2}{M_i} \\\\right ] \\\\pd{p_\\\\epsilon}\\\\nonumber\\\\\\\\\\n i L_{\\\\bo{p}} &=& \\\\sum_{i,\\\\alpha} \\\\bo{f}_{i,\\\\alpha} \\\\pd{\\\\bo{p}_{i,\\\\alpha}}\\\\nonumber\\\\\\\\\\n i L_{\\\\bo{q},MTK} &=& \\\\sum_{i,\\\\alpha} \\\\frac{p_\\\\epsilon}{W} \\\\bo{R}_i \\\\pd{\\\\bo{r}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\\\n i L_{\\\\bo{q}} &=& \\\\frac{\\\\bo{p}_{i,\\\\alpha}}{m_{i,\\\\alpha}} \\\\pd{\\\\bo{r}_{i,\\\\alpha}}\\\\nonumber\\\\\\\\\\n i L_{\\\\epsilon} &=& \\\\frac{p_\\\\epsilon}{W} 3 V \\\\pd{V}\\n\\\\end{eqnarray}\\n\\\\section{Integration of the equations of motion}\\nTo integrate the equations of motion we apply the following decomposition of the operator $e^{i \\\\hat{L}}$:\\n\\\\begin{eqnarray}\\n e^{i \\\\hat{L}} &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{p},MTK} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t}\\n e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}}\\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}}\\\\nonumber\\\\\\\\ \\n && e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n \\\\label{eq:trotdec} \\n\\\\end{eqnarray}\\nwhere all the operators are applied from left to right. Latter decomposition is based on the general approach which we briefly illustrate\\nin the following. If $iL_g$ a generic operator which can be split as follows:\\n\\\\begin{eqnarray}\\n i L_g = i L_1 + i L_2 + i L_3\\n\\\\end{eqnarray}\\none can use a first Trotter decomposition of the time evolution operator $e^{i L_g \\\\Delta t}$ as follows:\\n\\n\\\\begin{equation}\\n e^{iL_g \\\\Delta t } \\\\approx e^{i L_3 \\\\frac{\\\\Delta t}{2}} e^{(iL_1 + iL_2) \\\\Delta t} e^{i L_3 \\\\frac{\\\\Delta t}{2} } \\n\\\\end{equation}\\n\\nand one can employ a further Trotter decomposition for the operator \\n $e^{(iL_1 + iL_2) \\\\Delta t}$, i.e.\\n \\\\begin{eqnarray}\\n e^{(iL_1 + iL_2) \\\\Delta t} \\\\approx e^{iL_2 \\\\frac{\\\\Delta t}{2}} e^{i L_1 \\\\Delta t} e^{i L_2 \\\\dtd}\\n \\\\end{eqnarray}\\n\\nso that we have for the operator $e^{i L_g \\\\Delta t}$:\\n\\n\\\\begin{equation}\\n e^{iL_g \\\\Delta t} \\\\approx e^{iL_3 \\\\dtd} e^{iL_2 d\\\\dtd} e^{i L_1 \\\\Delta t} e^{i L_2 \\\\dtd} e^{i L_3 \\\\dtd} \\n \\\\label{eq:trotterMF}\\n\\\\end{equation}\\nThe operators $e^{i L_{NHP} \\\\Delta t/2}$ and $ e^{i L_{NHB} \\\\Delta t/2} $ can be factorized by leveraging the Trotter theorem, i.e.~if \\nwe define the following two operators:\\n\\\\begin{eqnarray}\\n S_{NHP}(\\\\delta) &=& e^{ \\\\frac{\\\\delta}{2} G_{\\\\eta,M_\\\\eta} \\\\pd{p_{\\\\eta_{M_\\\\eta}}} } \\n \\\\prod_{j=M_{\\\\eta}}^{1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{\\\\eta_j} \\\\pd{p_{\\\\eta_j}}} e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n \\\\right ) \\\\prod_{i=1}^{N_{at}} e^{\\\\delta \\\\frac{p_{\\\\eta_1}}{Q_{\\\\eta,1}} \\\\bo{p}_{i}\\\\cdot \\\\pd{\\\\bo{p}_i} }\\\\nonumber\\\\\\\\\\n && \\\\prod_{j=1}^{M_{\\\\eta}} e^{- \\\\delta\\\\frac{p_{\\\\eta_j}}{Q_{\\\\eta,j}} \\\\pd{\\\\eta_j} }\\n \\\\prod^{M_{\\\\eta}}_{j=1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{_j}\\\\pd{p_{\\\\eta_j}}} \\n e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\eta_{j+1}}}{Q_{\\\\eta,j+1}} \\\\pd{p_{\\\\eta_j}}} \\n \\\\right ) e^{ \\\\frac{\\\\delta}{2} G_{\\\\eta,M_\\\\eta} \\\\pd{p_{\\\\eta_{M_\\\\eta}}} }\\\\\\\\\\n S_{NHB}(\\\\delta) &=& e^{ \\\\frac{\\\\delta}{2} G_{,M_\\\\xi} \\\\pd{p_{_{M_\\\\xi}}} } \\n \\\\prod_{j=M_{\\\\xi}}^{1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{\\\\xi_j} \\\\pd{p_{\\\\xi_j}}} e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n \\\\right ) e^{\\\\delta \\\\frac{p_{\\\\xi_1}}{Q_{\\\\xi,1}} p_{\\\\epsilon} \\\\pd{p_\\\\epsilon} }\\\\nonumber\\\\\\\\\\n && \\\\prod_{j=1}^{M_{\\\\xi}} e^{-\\\\delta \\\\frac{p_{\\\\xi_j}}{Q_{\\\\xi,j}} \\\\pd{\\\\xi_j} }\\n \\\\prod^{M_{\\\\xi}}_{j=1} \\\\left ( e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n e^{ \\\\frac{\\\\delta}{2} G_{\\\\xi_j}\\\\pd{p_{\\\\xi_j}}} \\n e^{-\\\\frac{\\\\delta}{4} \\\\frac{p_{\\\\xi_{j+1}}}{Q_{\\\\xi,j+1}} \\\\pd{p_{\\\\xi_j}}} \\n \\\\right ) e^{ \\\\frac{\\\\delta}{2} G_{\\\\xi,{M_\\\\xi}} \\\\pd{p_{\\\\xi_{M_\\\\xi}}} }\\n \\\\label{eq:SNH}\\n\\\\end{eqnarray}\\none has:\\n\\\\begin{eqnarray}\\n e^{i L_{NHP} \\\\Delta t/2} &=& S_{NHP}(\\\\Delta t/2) + \\\\mathcal{O}({\\\\Delta t}^3) \\\\label{eq:trdecomNHP} \\\\\\\\ \\n e^{i L_{NHB} \\\\Delta t/2} &=& S_{NHB}(\\\\Delta t/2) + \\\\mathcal{O}({\\\\Delta t}^3)\\n \\\\label{eq:trdecomNHB}\\n\\\\end{eqnarray}\\n\\n$S_{NHP}$ and $S_{NHB}$ are called primitive decompositions of the operators $e^{i L_{NHP} \\\\Delta t/2}$ \\nand $ e^{i L_{NHB} \\\\Delta t/2} $ and they can be obtained by nothing that if one has a generic operator $i L_g$, \\nsuch that\\n\\\\begin{equation}\\n i \\\\wtL_g = i \\\\wtL_1 + i \\\\wtL_2 + i \\\\wtL_3 + i \\\\wtL_4\\n\\\\end{equation}\\nresorting to the same Trotter factorization used to arrive at Eq.~\\\\eqref{eq:trotterMF}\\nwith $iL_1 = i \\\\wtL$, $iL_2 = i \\\\wtL_2 + i\\\\wtL_3 $ and $ i L_3 = i \\\\wtL_4$, one obtains:\\n\\\\begin{eqnarray}\\n e^{i \\\\wtL_g \\\\Delta t} \\\\approx e^{i \\\\wtL_4 \\\\dtd} e^{ i (\\\\wtL_2 + i \\\\wtL_3) \\\\dtd} e^{ i \\\\wtL_1 \\\\Delta t}\\n e^{ i (\\\\wtL_2 + i \\\\wtL_3) \\\\dtd} e^{i \\\\wtL_4 \\\\dtd}\\n\\\\end{eqnarray}\\nand by factorizing again the operators $e^{ i (\\\\wtL_2 + i \\\\wtL_3) \\\\dtd}$, i.e.\\n\\\\begin{equation}\\n e^{(i\\\\wtL_2 + i \\\\wtL_3) \\\\dtd} \\\\approx e^{i \\\\wtL_3 \\\\dtq } e^{i \\\\wtL_2 \\\\dtd} e^{i \\\\wtL_3 \\\\dtq} \\n\\\\end{equation}\\none finally has:\\n\\\\begin{equation}\\n e^{i \\\\wtL_g \\\\Delta t} \\\\approx e^{i \\\\wtL_4 \\\\dtd} e^{i \\\\wtL_3 \\\\dtq } e^{i \\\\wtL_2 \\\\dtd} e^{i \\\\wtL_3 \\\\dtq} \\n e^{i \\\\wtL_4 \\\\dtd}\\n\\\\end{equation}\\n\\n% In Eqs.~\\\\eqref{eq:SYdecomNHP} and \\\\eqref{eq:SYdecomNHB} the weights $w_k$ can be chosen in such\\n% a way that an higher order integration scheme can be obtained. For example for a fourth-order scheme $n_{sy}=3$ and \\n% \\\\begin{eqnarray}\\n% w_1 &=& w_3 = \\\\frac{1}{2 - 2^{1/3}}\\\\nonumber\\\\\\\\\\n% w_2 &=& 1 - (w_1+w_3)\\n% \\\\label{eq:syweights}\\n% \\\\end{eqnarray}\\n% For a sixth-order scheme $n_{sy}=7$ and the weights are specified numerically~\\\\cite{tuckermanBook,Yoshida1990}. \\n% Note that if $n_{sy}=1$ and $w_1=1$ one recover the usual Trotter factorization.\\n\\nIn a similar fashion we can treat the operators $i L_{\\\\bo{q},MTK}$ and $i L_{\\\\bo{p},MTK}$. \\nIndeed, if we factorize the latter two operators as follows:\\n\\\\begin{eqnarray}\\n L_{\\\\bo{q},MTK} &=& \\\\sum_i L_{\\\\bo{q},MTK}^{(i)} = \\\\sum_i L_{\\\\bo{q},MTK0}^{(i)} + L_{\\\\bo{q},MTK1}^{(i)} \\\\nonumber\\\\\\\\ \\n L_{\\\\bo{p},MTK} &=& \\\\sum_i L_{\\\\bo{p},MTK}^{(i)} = \\\\sum_i L_{\\\\bo{p},MTK0}^{(i)} + L_{\\\\bo{p},MTK1}^{(i)} \\n \\\\label{eq:MTKop}\\n\\\\end{eqnarray}\\nwhere:\\n\\\\begin{eqnarray}\\n i L_{\\\\bo{q},MTK0}^{(i)} &=& \\\\sum_\\\\alpha \\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i}\\\\bo{r}_{i,\\\\alpha} \\\\pd{\\\\bo{r}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\ \\n i L_{\\\\bo{q},MTK1}^{(i)} &=& \\\\frac{p_\\\\epsilon}{W} \\\\sum_{\\\\ontop{\\\\alpha,\\\\alpha\'}{\\\\alpha\'\\\\neq\\\\alpha}}\\\\frac{m_{i,\\\\alpha\'}\\\\bo{r}_{i,\\\\alpha\'}}{M_i} \\\\pd{\\\\bo{r}_{i,\\\\alpha}} \\n \\\\label{eqMTKqsplit}\\n\\\\end{eqnarray}\\nand\\n\\\\begin{eqnarray}\\n i L_{\\\\bo{p},MTK0}^{(i)} &=& \\\\sum_{\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}}{M_i}\\\\bo{p}_{i,\\\\alpha} \\\\pd{\\\\bo{p}_{i,\\\\alpha}} \\\\nonumber\\\\\\\\ \\n i L_{\\\\bo{p},MTK1}^{(i)} &=& \\\\frac{p_\\\\epsilon}{W} \\\\sum_{\\\\ontop{\\\\alpha,\\\\alpha\'}{\\\\alpha\'\\\\neq\\\\alpha}}\\\\frac{m_{i,\\\\alpha}\\\\bo{p}_{i,\\\\alpha\'}}{M_i} \\\\pd{\\\\bo{p}_{i,\\\\alpha}} \\n \\\\label{eqMTKpsplit}\\n\\\\end{eqnarray}\\nagain we can employ the Trotter factorization to have:\\n\\\\begin{eqnarray}\\n e^{i L_{\\\\bo{q},MTK}^{(i)} \\\\Delta t/2} &=& S^{(i)}_{\\\\bo{q},MTK}(\\\\Delta t/2) + \\n \\\\mathcal{O}({\\\\Delta t}^3)\\\\label{eq:trdecomq} \\\\\\\\ \\n e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2} &=& S^{(i)}_{\\\\bo{p},MTK}(\\\\Delta t/2) + \\n \\\\mathcal{O}({\\\\Delta t}^3) \\n \\\\label{eq:trdecomp}\\n\\\\end{eqnarray}\\nwhere $S_{\\\\bo{q},MTK}(\\\\delta)$ and $S_{\\\\bo{p},MTK}(\\\\delta)$ are the following primitive factorizations of the operators \\n$e^{i L_{\\\\bo{q}}^{(i)}\\\\Delta t/2}$ and $e^{i L_{\\\\bo{p}}^{(i)} \\\\Delta t/2}$ respectively: \\n\\\\begin{eqnarray}\\n S^{(i)}_{\\\\bo{q},MTK}(\\\\delta) &=& \\\\prod_{\\\\alpha=1}^{n_i} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha\'}\\\\bo{r}_{i,\\\\alpha\'}}{M_i} \\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{r}_{i,\\\\alpha}}} \\n e^{i L_{\\\\bo{q},MTK0}^{(i)}\\\\delta } \\n \\\\prod_{\\\\alpha=n_i}^{1} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W}\\\\frac{m_{i,\\\\alpha\'}\\\\bo{r}_{i,\\\\alpha\'}}{M_i}\\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{r}_{i,\\\\alpha}}} \\n \\\\nonumber\\\\\\\\\\n S^{(i)}_{\\\\bo{p},MTK}(\\\\delta) &=& \\\\prod_{\\\\alpha=1}^{n_i} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}\\\\bo{p}_{i,\\\\alpha\'}}{M_i} \\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{p}_{i,\\\\alpha}}} \\n e^{i L_{\\\\bo{p},MTK0}^{(i)} \\\\delta} \\n \\\\prod_{\\\\alpha=n_i}^{1} e^{\\\\sum_{\\\\alpha\'\\\\neq\\\\alpha}\\\\frac{p_\\\\epsilon}{W} \\\\frac{m_{i,\\\\alpha}\\\\bo{p}_{i,\\\\alpha\'}}{M_i} \\\\frac{\\\\delta}{2} \\\\pd{\\\\bo{p}_{i,\\\\alpha}}}\\n \\\\label{eq:primfact} \\n\\\\end{eqnarray}\\n\\n\\\\noindent If the atoms of the $N$ molecules are subjected to the holonomic constraints:\\n\\\\begin{equation}\\n \\\\sigma_{i,\\\\beta}\\\\left (\\\\bo{r}_{i,1}\\\\ldots \\\\bo{r}_{i,n_i}\\\\right ) = 0\\n\\\\label{eq:constaints}\\n\\\\end{equation}\\nwhere $\\\\beta=1\\\\ldots l_i$ runs over all $l_i$ constraints of $i$-th molecule, the two operators $e^{i L_{\\\\bo{p}} \\\\Delta/2}$ \\nin Eq.~\\\\eqref{eq:trotdec} have to be modified to account for constraint forces\\nacting on each atom $(i,\\\\alpha)$ at time $t$ and $t+\\\\Delta t$. We assume that these forces can be expressed in terms of two different \\nsets of lagrangian multipliers ${\\\\lambda_{i,\\\\beta}(t)}$ and ${\\\\lambda_{i,\\\\beta}(t+\\\\Delta t)}$.\\nThe first operator from left in Eq.~\\\\eqref{eq:trotdec} evolves momenta according to forces evaluated at time $t$ and thus it becomes:\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{p}} \\\\Delta t/2} \\\\rightarrow e^{i L^{A}_{\\\\bo{p}} \\\\Delta t/2} = e^{\\\\left[\\\\sum_{i,\\\\alpha} \\\\bo{f}_{i,\\\\alpha} + \\n \\\\bo{f}^c_{i,\\\\alpha}(t) \\\\right] \\n \\\\pd{\\\\bo{p}_{i,\\\\alpha}}} \\n \\\\label{eq:op1constr}\\n\\\\end{equation}\\nwhere $\\\\bo{f}^c_{i,\\\\alpha}(t)$ is the sum of all constraint forces acting on the atom $(i,\\\\alpha)$, i.e. \\n\\\\begin{equation}\\n \\\\bo{f}^c_{i,\\\\alpha}(t) = \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t) \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{i,\\\\beta}(t)\\n\\\\end{equation}\\nand the second operator from left in Eq.~\\\\eqref{eq:trotdec} is modified as follows:\\n\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{p}} \\\\Delta t/2} \\\\rightarrow e^{i L^{B}_{\\\\bo{p}} \\\\Delta t/2} = e^{\\\\left[\\\\sum_{i,\\\\alpha} \\\\bo{f}_{i,\\\\alpha} \\n + \\\\bo{f}^c_{i,\\\\alpha}(t+\\\\Delta t) \\\\right]\\n \\\\pd{\\\\bo{p}_{i,\\\\alpha}}} \\n \\\\label{eq:op2constr}\\n\\\\end{equation}\\nwhere the constraint forces at time $t+\\\\Delta t$ are:\\n\\\\begin{equation}\\n \\\\bo{f}^c_{i,\\\\alpha}(t+\\\\Delta t) = \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t+\\\\Delta t) \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta (t+\\\\Delta t)\\n\\\\end{equation}\\nHence, the time evolution operator $e^{i \\\\hat{L}\\\\Delta t}$ becomes:\\n \\\\begin{eqnarray}\\n e^{i \\\\hat{L}} &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L^A_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t}\\n e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L^B_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}}\\\\nonumber\\\\\\\\ \\n && e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n \\\\label{eq:trotdec_constr}\\n\\\\end{eqnarray}\\n%We note that after the action of the operator $e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}}$ \\n%the correction of the position of an atom $(i,\\\\alpha)$ due to constraint forces \\n%depends only on constraint forces acting on $(i,\\\\alpha)$, i.e.\\\\@:\\n%\\\\begin{equation}\\n% e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\\\bo{r}_{i,\\\\alpha} = \\\\bo{r}_{i,\\\\alpha}^{nc} + \\\\frac{{\\\\Delta t}^2}{2} \\n% \\\\sum_\\\\beta \\\\tilde\\\\lambda_{i,\\\\beta}(t)\\n% \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{i,\\\\beta}(t)\\n%\\\\end{equation}\\n%where $\\\\bo{r}_{i,\\\\alpha}^{nc}$ is the position of the atom in absence of constraint forces and \\n%$\\\\tilde\\\\lambda_{i,\\\\beta}(t)$ are a set of unknown coefficients (see below on how to estimate them).\\n%Similarly, after the action of the operator $e^{i L_{\\\\bo{p},MTK} \\\\frac{\\\\Delta t}{2}}$ \\n%the correction of the velocity of an atom $(i,\\\\alpha)$ due to constraint forces \\n%depends only on constraint forces acting on $(i,\\\\alpha)$, i.e.\\\\@:\\n%\\\\begin{equation}\\n% e^{i L_{\\\\bo{p},MTK} \\\\frac{\\\\Delta t}{2}} \\\\bo{v}_{i,\\\\alpha} = \\\\bo{v}_{i,\\\\alpha}^{nc} + \\\\frac{\\\\Delta t}{2} \\\\sum_\\\\beta \\\\tilde\\\\lambda_{i,\\\\beta}(t+\\\\Delta t)\\n% \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{i,\\\\beta}(t)\\n%\\\\end{equation}\\n%where $\\\\bo{v}_{i,\\\\alpha}^{nc}$ is the velocity of the atom in absence of constraint forces.\\n\\nA convenient way to estimate the constraint forces is provided by \\nSHAKE algorithm~\\\\cite{Ciccotti1977,Ciccotti1982}. \\nWithin this approach, the full propagator is rewritten as follows:\\n\\\\begin{eqnarray}\\n e^{i \\\\hat{L}} &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}}\\n e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t} \\\\hat{\\\\mathcal{S}}_{r}\\n e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\\\nonumber\\\\\\\\\\n && \\\\hat{\\\\mathcal{S}}_{v}\\n e^{i L_{\\\\bo{p},MTK}\\\\frac{\\\\Delta t}{2}} e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}}\\n \\\\label{eq:trotdecshk}\\n\\\\end{eqnarray}\\nwhere the operator $\\\\hat{\\\\mathcal{S}}_{r}$ adjust positions and velocities of atoms to take into \\naccount the action of the constraint forces at time $t$ and the operator $\\\\hat{\\\\mathcal{S}}_{v}$ adjust the velocities \\naccording to constraint forces at time $t+\\\\Delta t$. Their action is described in detail in the following.\\n\\nThe operator $\\\\hat{\\\\mathcal{S}}_{r}$ adjusts first the positions of particles to fulfill the constraints in Eq.~\\\\eqref{eq:constaints}, then\\nupdates the velocities accordingly, i.e.~if the displacement of particle $(i,\\\\alpha)$ is $\\\\delta \\\\bo{r}_{i,\\\\alpha}$ the corresponding \\nvariation of its velocity is $\\\\delta \\\\bo{r}_{i,\\\\alpha}/\\\\Delta t$.\\nThe adjustment of the positions is achieved iteratively as in the standard SHAKE approach~\\\\cite{Ciccotti1977}, where\\nthe bonds, which always involve pair of atoms, are adjusted one by one until convergence is reached. \\n%We note that the displacement of two atoms involved in a bond\\n%has to occur along a direction calculated according to the positions of particles prior to the application of the operator\\n%$e^{i L_{\\\\bo{q}} \\\\Delta t}$. \\n%The calculation of the constraint forces $\\\\bo{f}^c_{i,\\\\alpha}(t)$ for updating the velocities \\n%requires some care.\\nIf we define the states $\\\\Gamma $, $\\\\Gamma\'$ and $\\\\Gamma\'\'$ as follows:\\n\\\\begin{eqnarray}\\n \\\\Gamma\' &=& e^{i L_{NHP} \\\\frac{\\\\Delta t}{2}} e^{i L_{NHB} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\n e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} \\\\Gamma(t) \\\\\\\\\\n \\\\Gamma\'\' &=& e^{i L_{\\\\bo{q}}\\\\Delta t} e^{i L_{\\\\epsilon}\\\\Delta t} \\\\, \\\\Gamma\'\\\\nonumber\\\\\\\\\\n %\\\\Gamma\'\' &=& e^{i L_{\\\\bo{q},MTK} \\\\Delta t /2}\\\\, \\\\Gamma\'\\n \\\\label{eq:rnoc}\\n\\\\end{eqnarray}\\nafter the application of $\\\\hat{\\\\mathcal{S}}_r$ the positions in state $\\\\Gamma\'\'$ will be adjusted\\nto fulfill all the holonomic constraints in Eq.~\\\\eqref{eq:constaints}, i.e.~if\\n\\\\begin{equation}\\n \\\\Gamma\'\'\' = \\\\hat{\\\\mathcal{S}}_{r} \\\\,\\\\Gamma\'\'% chktex 23 \\n \\\\label{eq:rp} \\n\\\\end{equation}\\nthe positions in state $\\\\Gamma\'\'\'$ fulfill the constraints in Eq.~\\\\eqref{eq:constaints}. % chktex 23\\n\\nLet\'s indicate with $\\\\bo{r}_{i,\\\\alpha}\'$, $\\\\bo{r}_{i,\\\\alpha}\'\'$ and $\\\\bo{r}_{i,\\\\alpha}\'\'\'$ \\nthe positions corresponding to states $\\\\Gamma\'$, $\\\\Gamma\'\'$ and $\\\\Gamma\'\'\'$ respectively.\\nThe operator $ \\\\hat{\\\\mathcal{S}}_{r} $ adjusts the positions by finding a set of lagrangian multipliers\\nsuch that:\\n\\\\begin{equation}\\n \\\\hat{\\\\mathcal{S}}_{r} \\\\bo{r}_{i,\\\\alpha}\'\' = \\\\bo{r}_{i,\\\\alpha}\'\' + \\\\frac{{\\\\Delta t}^2}{2} \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t) \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta (t)\\n \\\\label{eq:Sract}\\n\\\\end{equation}\\nwhere we note that $\\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta(t)$ must be calculated by using positions in state $\\\\Gamma\'$.\\nWe note that the operator $\\\\hat{S}_r$ can be safely applied before the operator $ e^{i L_{\\\\bo{q},MTK} \\\\frac{\\\\Delta t}{2}} $,\\nsince latter one does not alter the bonds, indeed:\\n\\\\begin{equation}\\n %\\\\dot{\\\\bo{r}}_{i,\\\\alpha} - \\\\dot{\\\\bo{r}}_{i,\\\\alpha\'} = \\n i L_{\\\\bo{q}}^{(i),MTK} (\\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'}) = 0\\n\\\\end{equation}\\nhence\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{q},MTK}^{(i)} \\\\frac{\\\\Delta t}{2}} (\\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'} ) = \\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'}\\n \\\\label{eq:distinv}\\n\\\\end{equation}\\nand from Eq.~\\\\eqref{eq:SYdecomq} one has that:\\n\\\\begin{equation}\\n \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{q},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n (\\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'} ) = \\n \\\\bo{r}\'\'\'_{i,\\\\alpha} - \\\\bo{r}\'\'\'_{i,\\\\alpha\'} + \\\\mathcal{O}({\\\\Delta t}^3)\\n \\\\label{eq:rconstrconserv}\\n\\\\end{equation}\\n%We note that the quantities $\\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_\\\\beta(t)$ are calculated by using positions \\n%in state $\\\\Gamma_0$.\\n\\nAs in the case of operator $\\\\hat{S}_r$ the adjustments of the velocities through the operator $\\\\hat{\\\\mathcal{S}}_{v}$\\ncan be viewed as the result of the action of constraint forces $\\\\bo{f}^c_{i,\\\\alpha}(t+\\\\Delta t)$ as the ones in\\nthe operator $e^{i L^B_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}}$ in Eq.~\\\\eqref{eq:trotdec_constr}.\\nIf we define the states $\\\\Gamma^{\\\\Romannum{4}}$ and $\\\\Gamma^{\\\\Romannum{5}}$ as follows:\\n\\\\begin{eqnarray}\\n \\\\Gamma^{\\\\Romannum{4}} &=& e^{i L_{\\\\bo{q}, MTK} \\\\frac{\\\\Delta t}{2}} \\\\Gamma\'\'\'\\\\\\\\\\n \\\\Gamma^{\\\\Romannum{5}} &=&e^{i L_{\\\\bo{p}} \\\\frac{\\\\Delta t}{2}} e^{i L_{p_{\\\\epsilon}} \\\\frac{\\\\Delta t}{2}} \\\\Gamma^{\\\\Romannum{4}} \\n \\\\label{eq:gamma45}\\n\\\\end{eqnarray}\\nthe action of the operator $\\\\hat{S}_v$ can be expressed as:\\n\\\\begin{equation}\\n \\\\Gamma^{\\\\Romannum{6}} = \\\\hat{S}_v \\\\Gamma^{\\\\Romannum{5}} = \\\\bo{v}^{\\\\Romannum{5}}_{i,\\\\alpha} + \\\\frac{{\\\\Delta t}^2}{2} \\\\sum_\\\\beta \\\\lambda_{i,\\\\beta}(t+\\\\Delta t)\\n \\\\nabla_{\\\\bo{r}_{i,\\\\alpha}} \\\\sigma_{\\\\beta}(t+\\\\Delta t)\\n\\\\end{equation}\\nwhere we note that $\\\\sigma_{\\\\beta}(t+\\\\Delta t)$ is calculated by using positions in state $\\\\Gamma^{\\\\Romannum{4}}$.\\nWe note that the operator $\\\\hat{S}_v$ can be safely applied after the operator $e^{i L_{\\\\bo{p}}^{(i)} \\\\Delta t/2}$, since latter\\none does not alter the relative velocity of atoms belonging to the same molecule, i.e.: \\n\\\\begin{equation}\\n i L_{\\\\bo{p},MTK}^{(i)} (\\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'}) = 0\\n\\\\end{equation}\\nso that\\n\\\\begin{equation}\\n e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\frac{\\\\Delta t}{2}} (\\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'} ) = \\n \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'}\\n \\\\label{eq:velinv}\\n\\\\end{equation}\\nand according to Eq.~\\\\eqref{eq:SYdecomp}:\\n\\\\begin{equation}\\n \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{p},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n (\\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\n \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'} ) = \\n \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha} - \\\\bo{v}^{\\\\Romannum{6}}_{i,\\\\alpha\'} + \\\\mathcal{O}({\\\\Delta t}^3)\\n \\\\label{eq:vconstrconserv}\\n\\\\end{equation}\\n\\n\\\\subsection{Suzuki-Yoshida decompositions}\\n\\nThe factorization of operators $e^{i L_{NHP}\\\\Delta t/2}$ and $e^{i L_{NHB}\\\\Delta t/2}$ in Eqs.~\\\\eqref{eq:trdecomNHP} and~\\\\eqref{eq:trdecomNHB}\\nand can be further improved by using a Suzuki-Yoshida (SY) decomposition and introducing RESPA \\nas discussed in~\\\\cite{tuckermanBook,Suzuki1991,Yoshida1990}.\\nSince $S_{NHP}$ and $S_{NHB}$ are primitive factorizations of the operators $e^{i L_{NHP}}$ and $e^{i L_{NHB}}$ we can \\nwrite~\\\\cite{tuckermanBook}:\\n%applying the SY decomposition to the operators $e^{i L_{NHP}\\\\Delta t/2}$ and $e^{i L_{NHB} \\\\Delta t/2}$ and by also introducing \\n%a RESPA with $n$ steps~\\\\cite{tuckermanBook} one has:\\n\\\\begin{eqnarray}\\ne^{i L_{NHP} \\\\Delta t/2} &\\\\approx& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S_{NHP}(w_k \\\\Delta t/2n)\\\\right]}^{n}\\\\label{eq:SYdecomNHP} \\\\\\\\ \\n e^{i L_{NHB} \\\\Delta t/2} &\\\\approx& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S_{NHB}(w_k \\\\Delta t/2n)\\\\right]}^{n} \\n \\\\label{eq:SYdecomNHB}\\n\\\\end{eqnarray}\\n\\nSimilarly, we can apply the Suzuki-Yoshida decomposition to the operators $e^{i L_{\\\\bo{q},MTK}^{(i)}\\\\Delta t/2}$ and \\n$e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2}$ and we can also introduce a RESPA with $n$ steps.\\nIn this case the primitive decompositions of the operators $e^{i L_{\\\\bo{q},MTK}^{(i)}\\\\Delta t/2}$ and \\n$e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2}$ are $S^{(i)}_{\\\\bo{q},MTK}$ and $S^{(i)}_{\\\\bo{p},MTK}$ so that:\\n\\n\\\\begin{eqnarray}\\n e^{i L_{\\\\bo{q},MTK}^{(i)} \\\\Delta t/2} &=& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{q},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n + \\n \\\\mathcal{O}({\\\\Delta t}^p) \\\\label{eq:SYdecomp} \\\\\\\\ \\n e^{i L_{\\\\bo{p},MTK}^{(i)} \\\\Delta t/2} &\\\\approx& \\\\prod_{k=1}^{n_{sy}} {\\\\left[ S^{(i)}_{\\\\bo{p},MTK}(w_k \\\\Delta t/2n)\\\\right]}^n \\n + \\\\mathcal{O}({\\\\Delta t}^p)\\n \\\\label{eq:SYdecomq}\\n\\\\end{eqnarray}\\nwhere $p$ is the order of the SY decomposition.\\nFinally, we note that thanks to SY decomposition the error in Eqs.~\\\\eqref{eq:rconstrconserv} and~\\\\eqref{eq:vconstrconserv} \\nbecomes $\\\\mathcal{O}({\\\\Delta t}^p)$, i.e.~it significantly reduces.\\n\\n\\\\section{Results and Discussion}\\nXXX\\n\\n\\\\section{Conclusions}\\n\\nYYY\\n\\n%\\\\subsection{Acknowledgments}\\n%bla, bla\\\\ldots\\n\\\\bibliography{isobaric}\\n\\\\end{document}\\n","uri":"file:///Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex","version":1}}}' 2021-05-13 11:00:18,241 - DEBUG - RX: Received message: b'{"jsonrpc":"2.0","method":"textDocument/publishDiagnostics","params":{"diagnostics":[],"uri":"file:///Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex"}}' 2021-05-13 11:00:18,242 - DEBUG - RX: Received message: b'{"jsonrpc":"2.0","method":"textDocument/publishDiagnostics","params":{"diagnostics":[],"uri":"file:///Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex"}}' 127.0.0.1 - - [13/May/2021 11:00:18] "POST /semantic_completion_available HTTP/1.1" 200 4 127.0.0.1 - - [13/May/2021 11:00:18] "POST /receive_messages HTTP/1.1" 200 171 127.0.0.1 - - [13/May/2021 11:00:20] "POST /debug_info HTTP/1.1" 200 618 2021-05-13 11:00:20,115 - DEBUG - RX: Received message: b'{"jsonrpc":"2.0","method":"textDocument/publishDiagnostics","params":{"diagnostics":[],"uri":"file:///Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex"}}' 127.0.0.1 - - [13/May/2021 11:00:20] "POST /receive_messages HTTP/1.1" 200 86 127.0.0.1 - - [13/May/2021 11:00:22] "POST /receive_messages HTTP/1.1" 200 5 2021-05-13 11:00:26,155 - INFO - No support for ExecuteCommand command in server for latex 2021-05-13 11:00:26,155 - INFO - No support for FixIt command in server for latex 2021-05-13 11:00:26,155 - INFO - Found definitionProvider support for command GoToDefinition in latex 2021-05-13 11:00:26,155 - INFO - Found definitionProvider support for command GoToDeclaration in latex 2021-05-13 11:00:26,155 - INFO - Found definitionProvider support for command GoTo in latex 2021-05-13 11:00:26,155 - INFO - No support for GoToType command in server for latex 2021-05-13 11:00:26,155 - INFO - No support for GoToImplementation command in server for latex 2021-05-13 11:00:26,155 - INFO - Found referencesProvider support for command GoToReferences in latex 2021-05-13 11:00:26,156 - INFO - Found renameProvider support for command RefactorRename in latex 2021-05-13 11:00:26,156 - INFO - Found documentFormattingProvider support for command Format in latex 2021-05-13 11:00:26,156 - INFO - Found workspaceSymbolProvider support for command GoToSymbol in latex 2021-05-13 11:00:26,156 - INFO - Found documentSymbolProvider support for command GoToDocumentOutline in latex 2021-05-13 11:00:26,156 - INFO - Always supporting StopServer for latex 2021-05-13 11:00:26,156 - INFO - Always supporting RestartServer for latex 2021-05-13 11:00:26,156 - INFO - Always supporting GetHover for latex 127.0.0.1 - - [13/May/2021 11:00:26] "POST /defined_subcommands HTTP/1.1" 200 148 2021-05-13 11:00:26,159 - INFO - No support for ExecuteCommand command in server for latex 2021-05-13 11:00:26,160 - INFO - No support for FixIt command in server for latex 2021-05-13 11:00:26,160 - INFO - Found definitionProvider support for command GoToDefinition in latex 2021-05-13 11:00:26,160 - INFO - Found definitionProvider support for command GoToDeclaration in latex 2021-05-13 11:00:26,160 - INFO - Found definitionProvider support for command GoTo in latex 2021-05-13 11:00:26,160 - INFO - No support for GoToType command in server for latex 2021-05-13 11:00:26,160 - INFO - No support for GoToImplementation command in server for latex 2021-05-13 11:00:26,160 - INFO - Found referencesProvider support for command GoToReferences in latex 2021-05-13 11:00:26,160 - INFO - Found renameProvider support for command RefactorRename in latex 2021-05-13 11:00:26,160 - INFO - Found documentFormattingProvider support for command Format in latex 2021-05-13 11:00:26,160 - INFO - Found workspaceSymbolProvider support for command GoToSymbol in latex 2021-05-13 11:00:26,160 - INFO - Found documentSymbolProvider support for command GoToDocumentOutline in latex 2021-05-13 11:00:26,160 - INFO - Always supporting StopServer for latex 2021-05-13 11:00:26,160 - INFO - Always supporting RestartServer for latex 2021-05-13 11:00:26,160 - INFO - Always supporting GetHover for latex 2021-05-13 11:00:26,161 - DEBUG - Refreshing file /Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex: State is Open/action None 2021-05-13 11:00:26,161 - DEBUG - TX: Sending message: b'Content-Length: 186\r\n\r\n{"id":2,"jsonrpc":"2.0","method":"textDocument/hover","params":{"position":{"character":0,"line":1},"textDocument":{"uri":"file:///Users/demichel/PAPERS/ISOBARICMD/LATEX/isobaric.tex"}}}' 2021-05-13 11:00:26,164 - DEBUG - RX: Received message: b'{"jsonrpc":"2.0","result":null,"id":2}' Traceback (most recent call last): File "/Users/demichel/.vim/plugged/YouCompleteMe/third_party/ycmd/third_party/bottle/bottle.py", line 868, in _handle return route.call(**args) File "/Users/demichel/.vim/plugged/YouCompleteMe/third_party/ycmd/third_party/bottle/bottle.py", line 1748, in wrapper rv = callback(*a, **ka) File "/Users/demichel/.vim/plugged/YouCompleteMe/third_party/ycmd/ycmd/watchdog_plugin.py", line 97, in wrapper return callback( *args, **kwargs ) File "/Users/demichel/.vim/plugged/YouCompleteMe/third_party/ycmd/ycmd/hmac_plugin.py", line 62, in wrapper body = callback( *args, **kwargs ) File "/Users/demichel/.vim/plugged/YouCompleteMe/third_party/ycmd/ycmd/handlers.py", line 91, in RunCompleterCommand return _JsonResponse( completer.OnUserCommand( ``` ------------------ file ```/var/folders/gc/2x7dz9ps7m38rc7ts03c30ch0000gn/T/ycmd_64159_stdout_eajav69t.log```: ``` serving on http://1.0.0.127.in-addr.arpa:64159 ``` ## OS version, distribution, etc. mac osx big sur 11.3.1 > Include system information here. ## Output of build/install commands > Include link to a [gist][] containing the invocation and entire output of > `install.py` if reporting an installation issue. [cont]: https://github.com/ycm-core/YouCompleteMe/blob/master/CONTRIBUTING.md [code]: https://github.com/ycm-core/YouCompleteMe/blob/master/CODE_OF_CONDUCT.md [readme]: https://github.com/ycm-core/YouCompleteMe/blob/master/README.md [faq]: https://github.com/ycm-core/YouCompleteMe/wiki/FAQ [search]: https://www.google.com/search?q=site%3Ahttps%3A%2F%2Fgithub.com%2Fycm-core%2FYouCompleteMe%2Fissues%20python%20mac [gist]: https://gist.github.com/

[puremourning]

Please reproduce without additional plugins (you included vimtex in your "minimal" reproduction).

[cridemichel]

yes indeed vimtex is not necessary to reproduce the problem, I have just edited the first post updating the minimal vimrc file

[puremourning]

to reproduce the problem you can open any latex file with the following minimal vimrc file:

From CONTRIBUTING.md:

Write a step-by-step procedure that when performed repeatedly reproduces your issue. If we can't reproduce the issue, then we can't fix it. It's that simple.

Explain what you expected to happen, and what actually happened. This helps us understand if it is a bug, or just a misunderstanding of the behavior.

Create a test case for your issue. This is critical. Don't talk about how "when I have X in my file" or similar, create a file with X in it and put the contents inside code blocks in your issue description. Try to make this test file as small as possible. Don't just paste a huge, 500 line source file you were editing and present that as a test. Minimize the file so that the problem is reproduced with the smallest possible amount of test data.

Can you include an exact test file and tell us exactly what tot type to reproduce the issue including what you expect to happen and what actually happens.

[cridemichel]

I have just edited my first post including also a minimal latex file to reproduce the issue best C.

[puremourning]

thanks, we'll take a look.

[cridemichel]

thank you

[bstaletic]

Here's what I see:

1. RefactorRename works.
2. None of the GoTo* commands ever return a location. 2.1. Well, GoToSymbol and GoToDocumentOutline do find lines 13 and 15.
3. GetHover works on arguments of \usepackage (as in \usepackage{argument}
4. GetHover receives null from texlab on argument of \begin and \end. Or anywhere else except 3.

The texlive log is empty. I don't see how this is YCM's fault. Perhaps the user is expected to have some documentation installed?

[cridemichel]

but as in my first post coc.nvim is able to provide documentation, why?

[puremourning]

This repo has nothing to do with coc.vim, so it's not really relevant.

[bstaletic]

@cridemichel You could have read my guess.

Perhaps the user is expected to have some documentation installed?

You can dive into coc's source and see what it is doing. Or you can ask on the texlive's support channels.

[cridemichel]

This repo has nothing to do with coc.vim, so it's not really relevant.

thank you very much for pointing it out but I was aware of that, I was just noticing that the coc.nvim is able to provide documentation from the textlab LSP, hence I thought that somehow YCM was failing on this, but if you say that YCM is working properly I believe you, thank you very much for your kindness and support

[bstaletic]

Doesn't look like our bug. YCM can't do anything with a null response. This question is really for the texlab maintainers.

[cridemichel]

I noticed that if I add set noshowmode in the minimal vimrc, which you can find in my first post, I switch to insert mode and I press "\" (backslash), I get the following error in the command line in red:

AssertionError:

which I presume it is triggered by YCM. Is this also related to texlab?

[bstaletic]

That's a bug in texlab. https://github.com/latex-lsp/texlab/issues/413

[cridemichel]

ok thank you again!

[oblitum]

@cridemichel

I was just noticing that the coc.nvim is able to provide documentation from the textlab LSP, hence I thought that somehow YCM was failing on this, but if you say that YCM is working properly I believe you

As pointed on latex-lsp/texlab#413. YCM is "working properly" by not having implemented an LSP feature which isn't mandatory, so technically it's at LSP server's fault to not comply with such constrained clients.

If you really wish full LSP support you should go for coc.nvim or nvimlsp anyways, as both started with the mindset of fully supporting LSP, instead of what happened/happens with YCM, which just added/adds LSP as an afterthought.