QM/MM documentation in user guide
[namd.git] / ug / ug_qmmm.tex
1 % QM/MM
2 % Documentation copied from https://www.ks.uiuc.edu/Research/qmmm/.
3
4 \section{Hybrid QM/MM Simulations}
5 \label{section:qmmm}
6
7 % This brief introduction aims at providing some basic context 
8 % to the following description of capabilities and commands 
9 % available in NAMD's QM/MM interface.
10
11 Even though molecular mechanics (MM) force-fields are based on
12 quantum mechanical calculations and experimental observations,
13 only quantum mechanics (QM) can give a complete and accurate
14 understanding of many biochemical processes, particularly those
15 involving chemical reactions or charge redistribution.
16 Nevertheless, even with the advanced hardware technology available today,
17 the computational cost of studying nanosecond-long dynamics of entire
18 systems relying solely on QM methodologies is usually prohibitive.
19 A common route to circumvent this cost barrier is to confine the
20 QM formalism to a sub-region of a system and to include the effects
21 of the surrounding system through MM simulations,
22 leading to hybrid QM/MM simulations~\cite{SENN2009}.
23
24 NAMD's comprehensive QM/MM suite~\cite{MELO2018} was developed to provide
25 easy setup, visualization and analysis of QM/MM simulations through
26 the graphical user interface VMD/QwikMD~\cite{RIBE2016},
27 and a broad range of QM methods through NAMD's new ``QMForces" module.
28 The QM/MM interface in NAMD supports the simulation of many
29 independent QM regions, and smooth integration with a vast collection
30 of enhanced sampling methods. In hybrid QM/MM simulations,
31 NAMD offloads part of its standard force and energy calculations
32 to a QM program, either through native interfaces to
33 MOPAC~\cite{STEW90,MAIA2012} or ORCA~\cite{NEES2012},
34 or through a flexible generic interface requiring a wrapper script,
35 where exemplary Python wrappers are provided for Gaussian,
36 TeraChem and Q-CHEM.
37 Multiple QM-MM coupling schemes are implemented, allowing for both
38 mechanically and electrostatically embedded QM regions to be used
39 (see description in Nature Methods~\cite{MELO2018}).
40 QM/MM simulations require the same input files used for classical MD,
41 with additional options in the configuration file.
42 QM and MM atoms covalently bound are usually treated by redistributing
43 the MM atom's charge over its nearest MM neighbors and by capping
44 the QM atom with a hydrogen atom,
45 as shown in Figure~\ref{fig:hybrid_qmmm} for a solvated
46 tri-alanine QM/MM calculation using the NAMD/ORCA interface.
47 Tests of the QM/MM interface for accuracy, stability and performance,
48 are provided as supporting information in Nature Methods~\cite{MELO2018}. 
49
50 If employing NAMD QM/MM please cite:
51 \begin{quote}
52 NAMD goes quantum: An integrative suite for hybrid simulations.
53 Melo*, M. C. R.; Bernardi*, R. C.; Rudack T.; Scheurer, M.; Riplinger, C.;
54 Phillips, J. C.; Maia, J. D. C.; Rocha, G. D.; Ribeiro, J. V.;
55 Stone, J. E.; Neese, F.; Schulten, K.; Luthey-Schulten, Z.;
56 Nature Methods, 2018 (doi:10.1038/nmeth.4638)
57 \end{quote}
58
59 \begin{figure}[tbp]
60 \centering
61 \includegraphics[width=6in]{figures/hybrid_qmmm_diagram.png}
62 \caption[Hybrid QM/MM NAMD]{%
63 Graphical representation of NAMD-ORCA interconnection.
64 Only the contribution of MM charges beyond rmax are
65 calculated by NAMD (via PME), with the direct electrostatic
66 calculation performed by ORCA.
67 The image assumes the charge shift redistribution scheme,
68 where the partial charge of the linking MM atom is shifted
69 to its nearest MM neighbors.
70 }
71 \label{fig:hybrid_qmmm}
72 \end{figure}
73
74
75 \subsection{Division of Labor}
76
77 The basic idea behind a hybrid QM/MM simulation in NAMD is to use 
78 a classical force field to treat the classical atoms in the system 
79 (or ``MM atoms"), and pass the information that describes the quantum 
80 atoms in the system (or ``QM atoms") to a Quantum Chemistry (QC) software, 
81 which is expected to produce gradients for all QM atoms, as well as 
82 the total energy of the QM region (and optionally partial charges). 
83 All bonded and non-bonded interactions among MM atoms are handled 
84 by NAMD's force field. Similarly, all bonded and non-bonded interactions 
85 among QM atoms are handled by the QC software in its chosen theory level. 
86 Treatment of covalent bonds between QM and MM atoms will be described 
87 in a following section.
88
89 The non-bonded interactions between QM and MM atoms are handled differently, 
90 and can be modified and regulated by the user. 
91 Van der Waals interactions are always calculated, and can be done 
92 using either the default force field parameters, or specific 
93 (user-defined) parameters for QM atoms. 
94 Parameter modifications for QM atoms have been proposed 
95 in order to compensate for over-polarization that these atoms 
96 may exhibit in hybrid QM/MM simulations. 
97 Larger van der Waals radii and/or shallower well depths 
98 should then be provided for all element types that occur 
99 among QM atoms (see the ``qmVdwParams" keyword).
100
101
102 \subsection{Mechanical and Electrostatic Embedding}
103
104 Electrostatic interactions between QM and MM atoms deserve 
105 a more detailed discussion due to the abundance and diversity 
106 of available alternatives. 
107 The first decision to be made is whether there will be electrostatic 
108 interactions between the two portions of a system, QM and MM. 
109 In the ``mechanical embedding" scheme, only positions and elements 
110 of atoms in the QM region are passed on to the chosen QC software 
111 for energy and force calculations. 
112 This way, QM and MM atoms share only van der  Waals interactions.
113
114 % Figure 1 for QM/MM
115 \begin{figure}[tbp]
116 \centering
117 \includegraphics[width=5in]{figures/OptionsDiagram.png}
118 \caption[Diagram of classical point charge options.]{%
119 Diagram of options that control the use and manipulation of 
120 classical point charges. Default values are indicated below their 
121 respective keyword. ``Cutoff" and ``SwitchDist" are keywords used in 
122 NAMD to configure the calculations of electrostatic and van der Waals 
123 interactions. 
124 }
125 \label{fig:qmmm_options}
126 \end{figure}
127
128 In the ``electrostatic embedding" scheme, on the other hand, 
129 the partial charges of MM atoms surrounding all QM atoms 
130 are used to approximate the electrostatic environment where QM atoms 
131 are found (the scheme is selected with the ``qmElecEmbed" keyword). 
132 See Figure~\ref{fig:qmmm_options}.
133 This process can be customized in a variety 
134 of ways, the first of which is deciding if a smoothing function will 
135 be used to avoid an abrupt decay in electrostatic force 
136 due to the cutoff used in the selection of surrounding point charges
137 (this option is activated with the ``qmSwitching" keyword).
138
139 Classical point charge utilization can be further customized 
140 by choosing which smoothing function will be used, 
141 and if the total charge of selected partial charges should be 
142 modified to (A) have a whole charge or (B) have a complementary 
143 charge to that of the QM region, so that the sum of charges 
144 from QM atoms and classical partial charges add to zero
145 (see Figure~\ref{fig:qmmm_options}).
146
147 With electrostatic embedding, QM atoms are influenced by the charges 
148 in the classical region. In order to balance the forces acting on the 
149 system, NAMD uses partial charges for the QM atoms to calculate 
150 the electrostatic interaction with classical point charges. 
151 There are two possibilities for the origin of the QM partial charges: 
152 the original partial charges found in the force field parameter files 
153 can be used, or updated partial charges can be gathered at each step 
154 from the QC software output 
155 (controllable through the ``qmChargeMode" keyword). 
156 The continuous update in charge distribution allows for a partial 
157 re-parameterization of the targeted molecule at each time step, 
158 which can lead to an improved description of the interactions 
159 of a ligand as it repositions over the surface of a protein, 
160 for example, or as it moves through a membrane.
161
162 In case PME is activated by the user, NAMD will automatically apply 
163 the necessary corrections to handle the QM region, allowing it to be 
164 influenced by long range interactions from the entire system.
165
166
167 \subsection{Covalent Bonds Divided by the QM/MM Barrier}
168
169 Hybrid QM/MM simulations of biomolecular systems often present 
170 situations where only a portion of a molecule should be treated 
171 quantum mechanically, usually to save computational resources 
172 since the cost of simulating QM regions rises rapidly with the 
173 number of simulated toms. In order to deal with chemical bonds 
174 that are split by the QM/MM division of the biomolecular system, 
175 that is, bonds that have one atom in the quantum (QM) region and 
176 another in the classical (MM) region (we will call these ``QM/MM bonds"), 
177 NAMD makes approximations to the molecular system in order to bridge 
178 differences in simulation type (QM vs.\ MM), and minimize errors involved 
179 in the QM/MM division of the system
180 (Figure~\ref{fig:qmmm_bond_treatment} A and B).
181
182 \begin{figure}[tbp]
183 \centering
184 \includegraphics[width=5in]{figures/QMMMBond_2bonds-01.png}
185 \caption[Treatment of QM/MM bonds]{%
186 Treatment of QM/MM bonds.
187 A) Illustration of all atoms in the vicinity of the QM/MM bond, 
188 colored by element: cyan for carbon, white for hydrogen, 
189 blue for nitrogen and red for oxygen. 
190 B) To the left, in blue, is the region that will be treated with 
191 the chosen QC software. To the right, in red, the region treated 
192 classically by NAMD. The bond in gray is the one crossing the 
193 QM/MM division. The atom marked as \textbf{QM1} is the quantum atom 
194 directly connected to the classical atom on the other side of the 
195 QM/MM division. Analogously, the atom marked as \textbf{MM1} is the 
196 classical atom directly connected to the quantum atom on the other 
197 side of the QM/MM division. Atoms marked as \textbf{MM2} are directly 
198 bonded to the \textbf{MM1} atom, and atoms marked \textbf{MM3} are 
199 directly bonded to \textbf{MM2} atoms. 
200 C) \textbf{Z1} method. Ignored partial charges are indicated in the image 
201 with a gray sphere taking the place of its respective classical atom. 
202 Directly adjacent to MM1 is a green sphere representing the link atom 
203 that is placed along the QM1-MM1 covalent bond. All remaining partial 
204 charges representing classical atoms that are passed on to the QC software
205 are indicated in purple spheres. 
206 D) \textbf{Z2} method.  
207 E) \textbf{Z3} method. 
208 F) RCD method. Virtual point charges, are represented in yellow spheres. 
209 The text indicates the total charge placed at each position, 
210 where \textbf{q} indicates the charge of the MM1 atom and \textbf{q2} 
211 represents the partial charge of the MM2 atom at that position. 
212 The yellow spheres at MM2 atom positions indicate their partial charge 
213 has been changed from its original value. 
214 G) CS method. Since in this case the virtual point charges are placed 
215 very close to the MM2 atom position, the yellow spheres representing 
216 them show significant overlapping.
217 }
218 \label{fig:qmmm_bond_treatment}
219 \end{figure}
220
221
222 \subsubsection{Link Atoms}
223
224 As previously mentioned, the information regarding atoms in the QM region 
225 is passed on to the chosen QC software, that is, their respective 
226 positions and element types, but in order to maintain (or approximate) 
227 the effect of the chemical bond between the QM atom and the MM atom, 
228 NAMD creates and places a link atom (usually a hydrogen) along the 
229 ``broken" QM/MM bond. The user can fine-tune this process by choosing 
230 the method of placement of the link atom and even the element of 
231 such atom (keywords ``qmBondDist" and ``qmLinkElement").
232
233 The introduction of the link atom will invariably place it very near 
234 the classical atom involved in the QM/MM bond, therefore the use 
235 and placement of partial charges from classical atoms becomes
236 highly relevant. Under the mechanical embedding scheme, the QC software 
237 only receives the atoms in the QM region and the link atoms created
238 to approximate QM/MM bonds, so no manipulation of partial charges is 
239 required. On the other hand, usual QM/MM simulations are done under 
240 the electrostatic embedding scheme, in which case the partial charges
241 of classical atoms involved in the QM/MM bonds and classical atoms 
242 directly connected to them require special treatment.
243
244
245 \subsubsection{Point Charge Alterations}
246
247 Several methods have been proposed to handle this situation, and the 
248 QM/MM interface developed here implements the most widely accepted ones.
249 One can be chosen using the ``qmBondScheme" keyword
250 (Figure~\ref{fig:qmmm_bond_treatment} C to G). 
251 In all implemented methods, the classical atom participating in the 
252 QM/MM bond (MM1 atom) does not have its partial charge passed on 
253 to the QC software, since this would create excessive repulsion 
254 (or attraction) on the link atom. This is, in fact, the entirety of 
255 the ``Z1" method: ignoring the partial charge of the MM1 atom. 
256 Analogously, ``Z2" and ``Z3" ignore all partial charges up to 
257 MM2 and MM3 atoms, respectively
258 (Figure~\ref{fig:qmmm_bond_treatment} C to E).
259
260 The Redistributed Charge and Dipole (RCD) method
261 (Figure~\ref{fig:qmmm_bond_treatment} F)
262 is more elaborate, as it rearranges the partial charge of the 
263 MM1 atom (indicated as $q$) so that the total charge of the region 
264 is maintained as well as the dipole moments of the bonds between 
265 MM1 and MM2 atoms. This is done by creating ``virtual" point charges, 
266 which are passed on to the QC software as if they represented 
267 partial charges of classical atoms. More specifically, the RCD method 
268 creates a virtual point charge in the middle of all MM1-MM2 bonds 
269 with a charge of $2q/n$, where $n$ is the number of 
270 MM2 atoms connected to MM1, and also subtracts a charge $q/n$ 
271 from each of the MM2 atoms, so that the total charge of the region 
272 remains constant while approximating the dipole moment of the 
273 MM1-MM2 bonds. This way there will be no point charge placed at the 
274 position of the MM1 atom, but its partial charge is not simply removed, 
275 it is redistributed.
276
277 A similar approach is taken by the Charge Shifting (CS) method 
278 (Figure~\ref{fig:qmmm_bond_treatment} G).
279 In this case, the MM1 partial charge is equally 
280 distributed across the MM2 atoms, and two virtual point charges 
281 are placed along the direction of the MM1-MM2 bond, one before 
282 the MM2 atom and one after, each one with a charge of $+q/n$ 
283 and $-q/n$ respectively. This method will also keep the 
284 total charge of the region constant while trying to preserve 
285 the local dipoles formed by all MM1-MM2 bonds.
286
287
288 \subsubsection{Link Atom Charge and Charge Groups}
289
290 Along with the gradient over all QM atoms, NAMD can also use 
291 the partial charge derived from the QC calculation to update 
292 the charge distribution of atoms in the QM region. 
293 When a QM/MM bond exists, however, part of the charge of the 
294 region will be placed on the link atom, and in order to keep the 
295 charge of the QM region constant, the link atom charge is 
296 re-distributed on the QM region. This seemingly simple mechanism 
297 can cause problems unless special care is be taken when deciding 
298 which bond will mark the division of QM and MM regions.
299
300 Many force fields divide the topologies of biomolecules in 
301 ``charge groups"
302 (Figure~\ref{fig:qmmm_chargegroups} A and B).
303 What this means is that 
304 not only will the partial charges of all atoms of a molecule 
305 add up to the whole number that represents the charge of the molecule, 
306 they will also add up to whole numbers in sub groups of atoms 
307 (look for the ``GROUP" statements in
308 \url{http://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/node24.html}
309 to see an example). 
310 Therefore, one needs to make sure that the chosen QM/MM bond(s) sits
311 in between ``charge groups", so the total sum of partial charges 
312 of atoms defining a QM region is a whole number. This is especially 
313 important in order to keep the total charge of the system constant. 
314 Since the QC calculation will always distribute a whole charge over 
315 all atoms of a system (QM atoms plus a link atom), if the partial charge 
316 of QM atoms is not initially a whole number, it will be forced into 
317 a whole number after the first QC step, where the charge of the link atom 
318 is distributed over the QM region. This will create a mismatch between 
319 QM and MM charges, changing the total charge of the entire system 
320 (QM plus MM regions).
321
322 \begin{figure}[tbp]
323 \centering
324 \includegraphics[width=5in]{figures/ChargeGroupDiagram.png}
325 \caption[Charge Groups and QM/MM Bonds]{%
326 Charge Groups and QM/MM Bonds.
327 A) Illustration of aspartate and the distribution of charge over its 
328 atoms as in CHARMM36 force field parameters. Circles in red indicate 
329 oxygen atoms, blue indicate nitrogen atoms, cyan for carbon atoms, 
330 and white for hydrogen atoms. ``Bond 1" indicates the peptide bond,
331 ``Bond 2" indicates the one between the alpha carbon and the peptide bond 
332 nitrogen, and ``Bond 3" the bond between the alpha carbon and the peptide 
333 bond carbon. 
334 B) Charge groups are indicated with dashed squares, 
335 along with their total charges. 
336 C) Depiction of the atoms in the QM region if Bonds 1 and 3 are used to 
337 separate it from the MM region. The total charge of QM region is $-1$. 
338 D) Depiction of QM region if the same is defined by Bonds 2 and 3. 
339 In this case, the total charge of QM region is $-1.16$.
340 }
341 \label{fig:qmmm_chargegroups}
342 \end{figure}
343
344 An example can be seen in Figure~\ref{fig:qmmm_chargegroups},
345 bonds 1 and 3 are chosen as the QM/MM bonds,
346 the charge distribution seen in Figure~\ref{fig:qmmm_chargegroups} C shows 
347 a whole charge for the QM region (and consequently for the MM region). 
348 Therefore, any charge placed on link atoms can be redistributed 
349 to the QM atoms with no change in total system charge. However, 
350 if bonds 2 and 3 are chosen for the QM/MM bond
351 (Figure~\ref{fig:qmmm_chargegroups} D), 
352 the charge of the MM region would be $+1.16$, while the charge 
353 of the QM region would be $-1.16$. Since the QC calculation would 
354 place a pre-determined whole charge on the region ($-1$, in this case), 
355 the updated charge of the QM region after the first simulation step 
356 would change the total charge of the system to $+0.16$, in this example.
357
358
359 \subsection{Custom Quantum Chemistry Software}
360
361 In order to offer the broad range of tools and technologies present 
362 in NAMD to all researchers who develop and/or employ specialized 
363 Quantum Chemistry tools, the QM/MM interface is prepared to utilize 
364 any QC tool that can be wrapped in a script that converts input and 
365 output files to specified formats. This flexible interface will 
366 improve development and testing of new tools, as well as their 
367 quick integration utilization in hybrid dynamics.
368
369 We currently provide 
370 in the \texttt{libs/qmmm/} directory
371 (and mirrored at
372 \url{http://www.ks.uiuc.edu/Research/qmmm/Scripts/})
373 Python wrapper scripts for GAUSSIAN, TeraChem, and Q-Chem.
374 Other wrapper scripts can be generated, based on these templates,
375 in any other language,
376 as long as they are provided to NAMD in an executable form.
377 Although natively supported,
378 we also provide a python wrapper script for ORCA,
379 with extended comments explaining the format in which NAMD will write 
380 data for the QC software and the format in which NAMD expects to find 
381 the results.
382
383
384 \subsection{Independent QM Regions}
385
386 Aiming at large macromolecular simulations that could take advantage 
387 of localized QM resolution, NAMD allows the user to set up multiple 
388 independent QM regions in the same molecular system. For example, 
389 one could study a multimeric complex that contains several active sites 
390 and have all active sites be calculated with a chosen QC software 
391 simultaneously (Figure~\ref{fig:qmmm_multiple_grid}).
392 Each active site would be calculated 
393 independently of all others, by its own execution of the QC software, 
394 keeping the calculation cost low and without impacting the overall 
395 efficiency of the simulation, since all QM regions would be calculated 
396 in parallel.
397
398 \begin{figure}[tbp]
399 \centering
400 \includegraphics[width=5in]{figures/MultipleQMDiagram.png}
401 \caption[Diagram of Multiple Grid Regions]{%
402 Diagram of Multiple QM Regions. 
403 The illustration depicts a hetero-hexameric complex 
404 (light and dark green triangles) that combine to create three active sites 
405 (in gray). Active sites are bound to a target molecule (purple triangle). 
406 The inner and outer dashed circles represent, respectively, the boundary 
407 of a QM region and the limits of the classical point charge shell around 
408 that region.
409 }
410 \label{fig:qmmm_multiple_grid}
411 \end{figure}
412
413 Identifying the different QM regions and which atoms belong to each 
414 one of them can be simply accomplished in the input PDB file, 
415 or in a dedicated PDB file (keyword ``qmParamPDB"). Since each region 
416 can contain different sets of molecules, their charges and multiplicities 
417 are indicated separately (see keywords ``qmCharge" and ``qmMult").
418
419 For simulations of large systems that are distributed across several 
420 computer nodes, one can control how many independent QM regions are 
421 calculated in each node. This would prevent large simulations from 
422 running out of memory if two or more large QM regions are placed 
423 in the same node (see keyword ``qmSimsPerNode").
424
425
426 \subsection{Keywords}
427
428 \begin{itemize}
429 \setlength{\itemsep}{0.4cm}
430
431 \item
432 \NAMDCONFWDEF{qmForces}{Calculate QM?}{yes or no}{no}{%
433 Turns on or off the QM calculations.
434 }
435
436 \item
437 \NAMDCONF{qmParamPDB}{Set QM atoms}{PDB file}{%
438 Name of an optional secondary PDB file where the OCCupancy
439 or BETA column has the indications for QM or MM atoms.
440 QM atoms should have an integer bigger than zero (0) and 
441 MM atoms should have zero as the beta or occupancy field. 
442 The same file may have indications for bonds between a QM atom 
443 and an MM atom. This should be provided in case the PDB file 
444 passed to the ``coordinates" keyword already has data on its 
445 beta or occupancy columns, such as when a SMD simulations 
446 is being performed.
447 }
448
449 \item
450 \NAMDCONF{qmColumn}{Which column?}{``beta" or ``occ"}{%
451 Indicates which column has the QM/MM field.
452 Required.
453 }
454
455 \item
456 \NAMDCONFWDEF{qmSimsPerNode}{Sims per node}{postive integer}{1}{%
457 Number of independent simultaneous QM simulations per node.
458 }
459
460 \item
461 \NAMDCONF{qmBondColumn}{Which bond column?}{``beta" or ``occ"}{%
462 Indicates which column has the QM/MM bond information.
463 This will tell NAMD which atoms are at the ends of a covalent bond
464 split by the QM/MM barrier, with one atom being quantum and
465 one being classical.
466 There is no default value.
467 If this parameter is provided,
468 NAMD will parse options for dealing with QM/MM bonds.
469 }
470
471 \item
472 \NAMDCONFWDEF{qmBondDist}{Use qmBondColumn value for distance?}%
473 {``on'' or ``off''}{off}{%
474 Indicates whether the value in the BondColumn will be used to define 
475 the distance between the QM atom and the Link Atom that will replace 
476 the MM atom in the QM system.
477 }
478
479 \item
480 \NAMDCONFWDEF{qmBondValueType}{Does qmBondColumn value give length or ratio?}%
481 {``len'' or ``ratio''}{len}{%
482 Indicates if the values in the BondColumn represent either 
483 the length (``len'') between the QM and link atoms or 
484 the ratio (``ratio'') between the QM-MM distance and 
485 the one which will be used as the QM-Link distance.
486 }
487
488 \item
489 \NAMDCONFWDEF{qmLinkElement}{Set link atom element}%
490 {string, for example ``4 9 Cl"}{H}{%
491 User defined link atom element. Two syntaxes are allowed:
492 if there is only one QM-MM bond, a string with the element symbol 
493 is allowed (such as ``H" or ``Cl"). If there are two or more bonds, 
494 the string needs to have the two atoms that compose the bond, 
495 and then the element (such as ``4 9 Cl"). 
496 The default element for all link atoms is hydrogen.
497 }
498
499 \item
500 \NAMDCONFWDEF{qmBondScheme}{Select QM-MM bond scheme}%
501 {``CS" or ``RCD" or ``Z1" or ``Z2" or ``Z3"}{CS}{%
502 Indicates what will be the treatment given to QM-MM bonds in terms 
503 of charge distribution and link atom creation and placement. 
504 CS: Charge Shift Scheme; RCD: Redistributed Charge and Dipole method; 
505 Z1: Only ignored MM1 partial charge, with no charge redistribution; 
506 Z2: Ignores MM1 and all MM2 partial charges, with no charge redistribution; 
507 Z3: Ignores MM1 and all MM2 and MM3 partial charges, 
508 with no charge redistribution.
509 }
510
511 \item
512 \NAMDCONFWDEF{qmElecEmbed}{Should point charges be used in QM?}%
513 {``on" or ``off"}{on}{%
514 Indicates if classical point charges should be used in QM calculations.
515 }
516
517 \item
518 \NAMDCONFWDEF{qmSwitching}{Use switching on point charges?}%
519 {``on" or ``off"}{off}{%
520 This will scale down the point charges representing the classical 
521 system as to replicate the switching procedure that NAMD applies 
522 to all charged interaction (see ``switching").
523 }
524
525 \item
526 \NAMDCONFWDEF{qmSwitchingType}{Set functional form of switching}%
527 {``switch" or ``shift"}{shift}{%
528 This option is used to decide which kind of function will be used 
529 to scale down point charges sent to QM calculations. 
530 SHIFT: This will "shift down" the entire shell of point charges 
531 so that electrostactic interactions reach zero at the cutoff distance. 
532 SWITCH: This will only change point charges in the sub-volume between 
533 the switchdist and cutoff distance, so that electrostactic interactions 
534 reach zero at the cutoff distance.
535 }
536
537 \item
538 \NAMDCONFWDEF{qmPointChargeScheme}{Set point charge scheme}%
539 {``none" or ``round" or ``zero"}{none}{%
540 This option allows the user to decide if and how the point charges 
541 presented to the QM system will be altered. NONE: Nothing will be done. 
542 ROUND: This will change the most distant point charges so that the 
543 total sum of point charges is a whole number. ZERO: This will adjust 
544 the most distant point charges so that the total sum of point charges 
545 is 0.
546 }
547
548 \item
549 \NAMDCONF{qmBaseDir}{Set directory for file I/O}{directory path}{%
550 This should be a fast read/write location, such as a RAM drive
551 (/dev/shm on most linux distros). The user needs to make sure this 
552 directory exists in the node(s) running the QM calculation(s).
553 }
554
555 \item
556 \NAMDCONFWDEF{qmReplaceAll}{Use only QM gradients for forces?}%
557 {``on" or ``off"}{off}{%
558 Indicates to NAMD that ALL forces form NAMD will be ignored and 
559 only the gradients from the QM software will be applied on the atoms. 
560 This IS NOT NECESSARY in any regular QM/MM simulation, and will prevent 
561 the use of any other feature from NAMD such as SMD.
562 }
563
564 \item
565 \NAMDCONFWDEF{qmVdwParams}{Modify type names for QM atoms?}%
566 {``on" or ``off"}{off}{%
567 The QM code will change all QM atoms' van der Waals types to "q"+element 
568 (e.g., all carbons will be qC and all hydrogens will be qH)
569 for VdW interactions. This means that a parameter file with 
570 epsilon and sigma values for these  atom types must be provided 
571 along with the regular parameter files. For example, if using 
572 CHARMM force field, the new file should be in CHARMM format.
573 }
574
575 \item
576 \NAMDCONFWDEF{qmPCStride}{Set stride for point charge}{integer}{1}{%
577 Sets a stride for new point charge determination. The same set of 
578 classical atoms will be sent to QM calculations as point charges, 
579 but with updated positions.
580 }
581
582 \item
583 \NAMDCONFWDEF{qmCustomPCSelection}{Provide custom point charge selection?}%
584 {``on" or ``off"}{off}{%
585 Indicates that one or more file(s) will be provided with a custom 
586 selection of point charges. Each file will have a selection for a single
587 QM group. This selection will be kept during the entire simulation.
588 }
589
590 \item
591 \NAMDCONF{qmCustomPCFile}{File for custom point charge selection}{PDB file}{%
592 The file will have, in the ``qmColumn", the same QM ID provided for
593 a single QM group. All other groups will have zero (0) in this column.
594 In the second column (beta or occupancy), the classical or quantum atoms
595 (from other QM regions) that need to be passed as point charges will be
596 identified by a non-zero number.
597
598 \vspace{-0.25em}
599 \textbf{Example/Format:}
600 \texttt{%
601 \vspace{-1em}
602 \small{
603 \begin{tabbing}
604 qmCustomPCFile \= system/system.customPC.1.pdb \\
605 qmCustomPCFile \> system/system.customPC.2.pdb \\
606 qmCustomPCFile \> system/system.customPC.3.pdb \\
607 qmCustomPCFile \> system/system.customPC.4.pdb
608 \end{tabbing}
609 } % \small
610 } % \texttt
611 } % \NAMDCONF
612
613 \item
614 \NAMDCONFWDEF{qmLiveSolventSel}{Keep track of solvent?}{``on" or ``off"}{off}{%
615 With Live Solvent Selection (LSS), NAMD will automatically keep track
616 of the solvent molecules for all QM Groups, and will exchange classical
617 solvent molecules with QM solvent molecules every "QMLSSFreq" steps.
618 }
619
620 \item
621 \NAMDCONFWDEF{qmLSSResname}{Set residue name for LSS}{residue name}{TIP3}{%
622 Indicates which residue name will be used in LSS.
623 }
624
625 \item
626 \NAMDCONFWDEF{qmLSSFreq}{Set frequency of LSS}%
627 {integer, multiple of stepspercycle}{100}{%
628 Frequency of LSS. Must be a multiple of stepspercycle.
629 }
630
631 \item
632 \NAMDCONFWDEF{qmLSSMode}{How solvent molecules are selected}%
633 {``dist" or ``COM"}{dist}{%
634 For LSS, this indicates how solvent molecules are selected.
635 In all cases, the closest solvent molecules are selected,
636 and if a classical solvent molecule is closer than a QM one,
637 they are swaped. DIST: This mode will use the smallest distance
638 between a solvent atom and a non-solvent QM atom to sort solvent
639 molecules. This is best used when the non-solvent QM atoms form
640 irregular volumes (when the COM is not very representatve),
641 and/or volumes with high solvent accessibility (such as a drug,
642 or a small peptide, in solution). COM: This mode will sort solvent
643 molecules based on Center Of Mass distance between the solvent COM
644 and the COM of a selection for each QM group (see ``qmLSSRef" keyword).
645 Best used with small QM regions that have limited solvent accessibility,
646 such as an active site.
647 }
648
649 \item
650 \NAMDCONF{qmLSSRef}{Which residues for COM of QM atoms?}{string}{%
651 This will indicate which residues are to be used in the determination
652 of the COM of non-solvent QM atoms. Only these atoms will be used to
653 determine the closest set of solvent molecules. The keyword takes a
654 string composed of the QM group ID, the segment name and the residue ID.
655
656 \vspace{-0.25em}
657 \textbf{Example/Format:}
658 \texttt{%
659 \vspace{-1em}
660 \small{
661 \begin{tabbing}
662 qmLSSRef \= "1 RP1 9" \\
663 qmLSSRef \> "2 RP1 3" \\
664 qmLSSRef \> "2 RP1 2" \\
665 qmLSSRef \> "3 AP1 9" \\
666 qmLSSRef \> "3 AP1 3" \\
667 qmLSSRef \> "4 AP1 9"
668 \end{tabbing}
669 } % \small
670 } % \texttt
671 } % \NAMDCONF
672
673 \item
674 \NAMDCONF{qmConfigLine}{Pass string to QM configuration}{string}{%
675 The string passed to "qmConfigLine" will be copied and pasted at the
676 very begining of the configuration file for the chosen QM software
677 if either ORCA or MOPAC are selected.
678
679 \vspace{-0.25em}
680 \textbf{Example/Format (QM/MM NAMD-ORCA):}
681 \texttt{%
682 \vspace{-1em}
683 \small{
684 \begin{tabbing}
685 qmConfigLine \= "! PM3 ENGRAD" \\
686 qmConfigLine \> "\%\%output PrintLevel Mini Print\textbackslash{[}P\_Mulliken\textbackslash{]} 1 Print\textbackslash{[}P\_AtCharges\_M\textbackslash{]} 1 end"
687 \end{tabbing}
688 } % \small
689 } % \texttt
690
691 \vspace{-1em}
692 \textbf{Example/Format (QM/MM NAMD-MOPAC):}
693 \texttt{%
694 \vspace{-1em}
695 \small{
696 \begin{tabbing}
697 qmConfigLine \= "PM7 XYZ T=2M 1SCF MOZYME CUTOFF=9.0 AUX LET GRAD QMMM GEO-OK" \\
698 qmConfigLine \> "Test System"
699 \end{tabbing}
700 } % \small
701 } % \texttt
702
703 } % \NAMDCONF
704
705 \item
706 \NAMDCONF{qmMult}{Set multiplicity of QM region}{string}{%
707 Multiplicity of the QM region. This is needed for proper construction
708 of ORCA's input file. Each string must be composed of the QM region ID
709 and its multiplicity.
710
711 \vspace{-0.25em}
712 \textbf{Example/Format:}
713 \texttt{%
714 \vspace{-1em}
715 \small{
716 \begin{tabbing}
717 qmMult \= "1 1" \\
718 qmMult \> "2 1" \\
719 qmMult \> "3 1" \\
720 qmMult \> "4 1"
721 \end{tabbing}
722 } % \small
723 } % \texttt
724 } % \NAMDCONF
725
726 \item
727 \NAMDCONF{qmCharge}{Set charge of each QM region}{string}{%
728 Indicates the charge of each QM region. If no charge is provided
729 for a QM region, NAMD calculates the total charge automatically
730 based on the given parameter set. Each string must be composed
731 of the QM region ID and its total charge.
732
733 \vspace{-0.25em}
734 \textbf{Example/Format:}
735 \texttt{%
736 \vspace{-1em}
737 \small{
738 \begin{tabbing}
739 qmCharge \= "1 1" \\
740 qmCharge \> "2 -1" \\
741 qmCharge \> "3 1" \\
742 qmCharge \> "4 -1"
743 \end{tabbing}
744 } % \small
745 } % \texttt
746 } % \NAMDCONF
747
748 \item
749 \NAMDCONF{qmSoftware}{Which QM software?}%
750 {``mopac" or ``orca" or ``custom"}{%
751 Required for QM/MM, this
752 indicates which QM software should be used. In case the user wants
753 to apply another QM code, this can be done by using the "custom"
754 qmSoftware. In this case, NAMD will call the executable defined in
755 the qmExecPath variable and will give it one argument: the full path
756 to the input file. 
757
758 INPUT: This input file will contain on the first line the number of
759 QM atoms (X) and the number of point charges in the file (Y, which
760 may be 0 or more), separated  by a space. The following X+Y lines
761 will have four (4) fields: X, Y and Z coordinates, and a fourth field
762 which will depend on the type of entry. For QM atoms, the  field will
763 contain the element of the QM atom. For point charge lines, the field
764 will contain the charge of the point charge.
765
766 OUTPUT: The expected output file whould be placed in the same directory
767 as the input file, and should be named \texttt{"*inputfile*.result"}
768 (meaning it will have the same path and name of the input file, plus the
769 suffix \texttt{".result"}). This file should have, on its first line,
770 the energy of the system, and on the following X lines (where X is the
771 number of QM atoms passed in the input file), four (4) fields: the x, y,
772 and z components of the TOTAL FORCE applied on that atom, and on the
773 fourth field, the charge of the atom. If the user indicates that charges
774 from the QM software should not be used (see ``qmChargeMode"), the fourth
775 field should have zeroes, but should not be empty. Energy should be
776 in Kcal/mol and forces in Kcal/mol/Angstrom.
777 }
778
779 \item
780 \NAMDCONF{qmExecPath}{Set path to QM executable}{path}{%
781 Required for QM/MM, this
782 indicates the path to the QM code executable.
783 }
784
785 \item
786 \NAMDCONF{qmSecProc}{Set path to secondary executable}{path}{%
787 Indicates a secondary executable that NAMD will call AFTER each
788 QM software execution for each QM group. The executable is called
789 with two arguments: the complete path and name of the input file
790 used for each QM software execution; and the simulation step. This
791 option can be used for an extra-processing at every step, e.g.,
792 for saving all QM outputs every step.
793 }
794
795 \item
796 \NAMDCONF{qmPrepProc}{Set path to initial executable}{path}{%
797 Indicates an executable that must be called BEFORE the FIRST QM
798 software execution of each QM group. The executable is called with
799 one argument: the complete path and name of the input file used
800 for each QM software execution. This can be used to setup a charge
801 distribution for a molecule of interest, for example.
802 }
803
804 \item
805 \NAMDCONFWDEF{qmChargeMode}{Set charge calculation mode}%
806 {``none" or ``mulliken" or ``chelpg"}{mulliken}{%
807 Charge calculation mode expected from the QM software. This
808 indicates if charges should be read from the QM software and
809 updated at every step, or if the original force field atom
810 charges should be used. In case you are using ORCA, two charge
811 options are allowed, Mulliken or CHELPG. We must know the kind
812 of charge requested by the user so that the proper format is expected,
813 read and processed. NONE: No charges are read from the QM software
814 output and the original force field charges are preserved.
815 MULLIKEN: This is the only other option for MOPAC and one possibility
816 for ORCA. In case you are using the custom QM software interface,
817 choose this option in order to use the charges calculated in the
818 custom QM software, irrespective of the actual theory used for
819 that calculation. CHELPG: This is a second possibility for ORCA.
820 }
821
822 \item
823 \NAMDCONFWDEF{qmOutStride}{Set frequency of QM charge output}{integer}%
824 {0 (not saving)}{%
825 Frequency of QM charge output. A dedicated DCD file will be created
826 to store the charge data for all QM atoms in the system. This
827 independent frequency allows the user to store whole-system data
828 at a larger stride to save time and space.
829 }
830
831 \item
832 \NAMDCONFWDEF{qmPositionOutStride}{Set frequency of QM-only position output}%
833 {integer}{0 (not saving)}{%
834 Frequency of QM-only position output. A dedicated DCD file will be
835 created to store the position data for all QM atoms in the system.
836 This independent frequency allows the user to store whole-system
837 data at a larger stride to save time and space.
838 }
839
840 \end{itemize}
841