QM/MM documentation in user guide 58/4758/4
authorDavid Hardy <dhardy@ks.uiuc.edu>
Wed, 31 Oct 2018 12:59:34 +0000 (07:59 -0500)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Thu, 8 Nov 2018 20:43:12 +0000 (14:43 -0600)
Documentation has been created from the research web page and paper.
Included introductory material with figures from website, pointer
to Python scripts, keyword definitions, and references.  Tested
building both the PDF and HTML.

Also have included Python and other relevant scripts into a new
lib/qmmm/ directory.

Change-Id: I0e9376500425af20b6d7d6bb6529b3728cd9897f

16 files changed:
lib/qmmm/cleanFolder.sh [new file with mode: 0755]
lib/qmmm/orcaviewer.tcl [new file with mode: 0644]
lib/qmmm/runORCA.py [new file with mode: 0755]
lib/qmmm/run_gaussian.py [new file with mode: 0755]
lib/qmmm/run_qchem.py [new file with mode: 0755]
lib/qmmm/run_terachem.py [new file with mode: 0755]
lib/qmmm/saveDipole.sh [new file with mode: 0755]
ug/Makefile
ug/figures/ChargeGroupDiagram.png [new file with mode: 0644]
ug/figures/MultipleQMDiagram.png [new file with mode: 0644]
ug/figures/OptionsDiagram.png [new file with mode: 0644]
ug/figures/QMMMBond_2bonds-01.png [new file with mode: 0644]
ug/figures/hybrid_qmmm_diagram.png [new file with mode: 0644]
ug/ug.bib
ug/ug.tex
ug/ug_qmmm.tex [new file with mode: 0644]

diff --git a/lib/qmmm/cleanFolder.sh b/lib/qmmm/cleanFolder.sh
new file mode 100755 (executable)
index 0000000..816ba88
--- /dev/null
@@ -0,0 +1,29 @@
+#!/bin/bash
+
+# If more than one independent QM region exists in the system, but only one QM
+# simulation per node is requested in the configuration file, all independent QM
+# regions will be ran in the same folder (the idea behind limiting the number of QM
+# regions being calculated per node is exactly saving RAM space). Some softwares,
+# however, may leave behind temporary files which are QM region-specific, which 
+# will lead to geometry errors when running multiple QM regions with one simulation
+# per node.
+#
+# This script erases temporary files left behing by ORCA so that two *different* and 
+# independent QM regions can be calculated in the same node (or single computer) 
+# in *sequence*, and NOT in *parallel*. Without the script, temporary files left 
+# behind after the first QM region was processed would be read by ORCA when 
+# calculating the second QM region, but since they are different, ORCA would crash
+# correctly claiming that the geometry in temporary files does not match the 
+# geometry in the input file.
+# 
+# This will NOT be more efficient, as independent QM regions could be calculated in
+# parallel, speeding up the overall time per simulation step, but the script serves 
+# as an example of what could be done with a Secondary Process.
+
+targDir=`dirname "$1"`
+
+echo "Input file: $1"
+echo "Target Dir: $targDir"
+
+rm $targDir/*.gbw
+rm $targDir/*.prop
\ No newline at end of file
diff --git a/lib/qmmm/orcaviewer.tcl b/lib/qmmm/orcaviewer.tcl
new file mode 100644 (file)
index 0000000..22f83b1
--- /dev/null
@@ -0,0 +1,282 @@
+package provide orcaviewer 1.0
+package require tablelist
+
+proc enabletrace {} {
+  global vmd_frame;
+  trace variable vmd_frame([molinfo top]) w updateFrame
+}
+
+proc disabletrace {} {
+  global vmd_frame;
+  trace vdelete vmd_frame([molinfo top]) w updateFrame
+}
+
+
+namespace eval OrcaViewer:: {
+  namespace export orcaviewer
+
+  variable w                                          ;# handle to main window
+  variable orcaFile   ""                  ;# output file
+  variable fileLoaded 0
+  variable orcaMol ""
+  variable orbitals {}
+  variable orbitalsInTable 0
+  variable lowestOrbitalIndex 0
+  variable reps {}
+  variable orbSpan 10
+}
+
+#
+# Create the window
+#
+proc OrcaViewer::orcaviewer {} {
+  variable w
+  variable orcaFile
+  variable fileLoaded
+  variable ::OrcaViewer::orbSpan
+
+  # If already initialized, just turn on
+  if { [winfo exists .orcaviewer] } {
+    wm deiconify $w
+    return
+  }
+
+  set w [toplevel ".orcaviewer"]
+  wm title $w "OrcaViewer"
+  wm resizable $w 0 0
+
+  frame $w.menubar -relief raised -bd 2 ;# frame for menubar
+  pack $w.menubar -padx 1 -fill x
+
+  menubutton $w.menubar.help -text Help -underline 0 -menu $w.menubar.help.menu
+  menubutton $w.menubar.file -text File -underline 0 -menu $w.menubar.file.menu
+
+  ## help menu
+  # menu $w.menubar.help.menu -tearoff no
+  # $w.menubar.help.menu add command -label "Help..." -command "vmd_open_url [string trimright [vmdinfo www] /]/plugins/textview"
+  # XXX - set menubutton width to avoid truncation in OS X
+  # $w.menubar.help config -width 5
+
+  menu $w.menubar.file.menu -tearoff no
+  $w.menubar.file.menu add command -label "Load" -command  OrcaViewer::loadFile
+  $w.menubar.file.menu add command -label "Clean" -command OrcaViewer::cleanUp
+  # $w.menubar.file.menu add command -label "Save" -command  TextView::savefile
+  # $w.menubar.file.menu add command -label "Save As" -command  TextView::saveasfile
+  $w.menubar.file config -width 5
+  pack $w.menubar.file -side left
+  pack $w.menubar.help -side right
+
+
+  ##
+  ## main window area
+  ##
+  wm resizable $w 1 1
+  wm geometry $w 500x200
+  # frame $w.main
+  scrollbar $w.scr -orient vertical -command [list $w.main  yview]
+  tablelist::tablelist $w.main -columns {0 "Orb. Number" 0 "Descr." 0 "Energy"} -stretch all -background white \
+    -selectmode extended \
+    -exportselection true \
+    -state normal -selectmode extended \
+    -labelrelief groove \
+    -labelbd 1 \
+    -yscrollcommand [list $w.scr set] \
+    -editstartcommand OrcaViewer::startOrbitalClick -editendcommand OrcaViewer::updateOrbitals
+
+    bind $w.main <<TablelistSelect>> {
+        set col [tablelist::getTablelistColumn %W]
+        set row [%W curselection]
+        puts $row
+        OrcaViewer::updateOrbitals $row $col
+    }
+  pack $w.main -fill both -expand 1 -side top
+  #$w.main insert end [list "first row" "another value"]
+  #$w.main insert end [list "another row" "bla bla"]
+
+  # label $w.txt.label -width 80 -relief sunken -bg White -textvariable TextView::textfile
+  # text $w.txt.text -bg White -bd 2 -yscrollcommand "$::TextView::w.txt.vscr set"
+  # scrollbar $w.txt.vscr -command "$::TextView::w.txt.text yview"
+  # pack $w.txt.label
+  # pack $w.txt.text $w.txt.vscr -side left -fill y
+  pack $w.menubar $w.main -side top -expand 1
+  entry $w.span -textvariable ::OrcaViewer::orbSpan -width 10
+  bind $w.span <Return> { OrcaViewer::fillOrblist }
+  label $w.spanLabel -text "Orbitals around HOMO: "
+  pack $w.main $w.span -expand 1 -side bottom
+  pack $w.spanLabel $w.span -expand 1 -side right
+}
+
+proc OrcaViewer::updateOrbitals {row col} \
+{
+  variable orbitals
+  variable orcaMol
+  variable lowestOrbitalIndex
+  variable reps
+  set orb [expr $row+$lowestOrbitalIndex+1]
+
+  if {[llength $orbitals] && ![llength $reps]} {
+    #puts "update orbital $row"
+    puts [lindex $orbitals [expr $row]]
+    mol color ColorID 0
+    mol representation Orbital 0.050000 $orb 0 0 0.125 1 0 0 0 1
+    mol selection all
+    mol material Glossy
+    mol addrep $orcaMol
+    mol color ColorID 3
+    mol representation Orbital -0.050000 $orb 0 0 0.125 1 0 0 0 1
+    mol addrep $orcaMol
+    set nr [molinfo $orcaMol get numreps]
+    set reps [list [mol repname $orcaMol [expr $nr - 1]] [mol repname $orcaMol [expr $nr - 2]]]
+  } elseif {[llength $orbitals] && [llength $reps]} {
+    set idx1 [mol repindex $orcaMol [lindex $reps 0]]
+    set idx2 [mol repindex $orcaMol [lindex $reps 1]]
+    mol modstyle $idx1 $orcaMol Orbital 0.050000 $orb 0 0 0.125 1 0 0 0 1
+    mol modstyle $idx2 $orcaMol Orbital -0.050000 $orb 0 0 0.125 1 0 0 0 1
+  }
+
+}
+
+proc OrcaViewer::getOrbitalDescr {index hindex} \
+{
+  set diff [expr $hindex-$index]
+  if {$diff == 0} {
+      return "HOMO"
+  } elseif {$diff > 0} {
+      return "HOMO-$diff"
+  } elseif {$diff < 0} {
+    if {$diff == -1} {
+      return "LUMO"
+    } else {
+      return "LUMO+[expr -$diff-1]"
+    }
+  }
+
+}
+
+proc OrcaViewer::assignPCtoUserField { method } {
+  variable orcaMol
+  set nf [molinfo $orcaMol get numframes]
+  set sel [atomselect $orcaMol all]
+  for {set i 0} {$i < $nf} {incr i} {
+    $sel frame $i
+    animate goto $i
+    set qmcharges [molinfo $orcaMol get qmcharges]
+    foreach var $qmcharges {
+      if {[lindex $var 0 0] == $method} {
+        $sel set user [lindex $var 0 1]
+        break
+      }
+    }
+  }
+  $sel delete
+}
+
+proc updateFrame { name element op } {
+  OrcaViewer::fillOrblist
+}
+
+proc OrcaViewer::fillOrblist { } \
+{
+  variable w
+  variable orcaMol
+  variable orbitals
+  variable orbitalsInTable
+  variable lowestOrbitalIndex
+  variable ::OrcaViewer::orbSpan
+  global vmd_frame
+
+  #puts $orbitalsInTable
+  $w.main delete 0 $orbitalsInTable
+  set orbitals {}
+
+  set energies [lindex [molinfo $orcaMol get orbenergies] 0 0]
+  set norbs [llength $energies]
+  #puts "norbs: $norbs"
+  set homo [molinfo $orcaMol get homo]
+
+  set orbitalsInTable 0
+  set lowestOrbitalIndex [expr $homo - $::OrcaViewer::orbSpan]
+  if {$lowestOrbitalIndex < 0} {
+    set lowestOrbitalIndex 0
+  }
+  set limit [expr $homo + $::OrcaViewer::orbSpan]
+  if {$limit > $norbs} {
+    set limit [expr $norbs - 1]
+  }
+  for {set i $lowestOrbitalIndex} {$i <= $limit} {incr i} {
+    lappend orbitals [list $i [OrcaViewer::getOrbitalDescr $i $homo] [lindex $energies $i]]
+    incr orbitalsInTable
+  }
+
+  foreach orbital $orbitals {
+    $w.main insert end $orbital
+  }
+
+}
+
+
+proc OrcaViewer::loadFile { } {
+  variable w
+  variable orcaFile
+  variable fileLoaded
+  variable orcaMol
+
+  set file_types {
+    {"All Files" * }
+  }
+
+  set orcaFile [tk_getOpenFile -filetypes $file_types \
+                -initialdir pwd \
+                -defaultextension .txt]
+
+  set rc [ catch { set fd [open $orcaFile "r"] } ]
+  if { $rc == 1} {
+    return
+  }
+
+  set fileLoaded 1
+  puts $orcaFile
+  puts $fileLoaded
+
+  set orcaMol [mol load orca $orcaFile]
+  puts $orcaMol
+  #OrcaViewer::fillOrblist
+
+  enabletrace
+  OrcaViewer::assignPCtoUserField "Mulliken"
+  color scale method BWR
+  animate goto 0
+
+  close $fd
+}
+
+proc OrcaViewer::cleanUp { } \
+{
+  variable w
+  variable orcaFile
+  variable fileLoaded
+  variable orcaMol
+  variable orbitalsInTable
+  variable reps
+  variable orbitals
+  variable lowestOrbitalIndex
+
+  mol delete $orcaMol
+  $w.main delete 0 $orbitalsInTable
+  set fileLoaded 0
+  set orcaMol ""
+  set orcaFile ""
+  set reps {}
+  set orbitals {}
+  set orbitalsInTable 0
+  set lowestOrbitalIndex 0
+}
+
+proc orcaviewer_tk {} {
+  OrcaViewer::orcaviewer
+  return $OrcaViewer::w
+}
+
+vmd_install_extension orcaviewer orcaviewer_tk "Visualization/Orca"
+orcaviewer_tk
+enabletrace
diff --git a/lib/qmmm/runORCA.py b/lib/qmmm/runORCA.py
new file mode 100755 (executable)
index 0000000..cf9d0f7
--- /dev/null
@@ -0,0 +1,344 @@
+#!/usr/bin/python3
+
+# by Marcelo Melo (melomcr@gmail.com)
+
+# If you would like to run a Quantum Chemistry software which is not supported by NAMD,
+# an interface script can be created in order to convert NAMD's input into the format
+# expected by the QC software, and subsequently convert the software's output inro the
+# format expected by NAMD.
+# 
+# This script creates an interface for ORCA, so that it can be ran using the custom 
+# QM software option provided by NAMD. It is an example script, since ORCA is supported
+# natively by NAMD's QM/MM interface, and is intended to be used by those wishing to 
+# test theyr own software in hybrid MD simulations through NAMD.
+
+#######  NAMD Files Format
+# INPUT: This input file will contain on the first line the number of QM atoms (X)
+# and the number of point charges in the file (Y, which may be 0 or more), separated 
+# by a space. The following X+Y lines will have four (4) fields: X, Y and Z coordinates,
+# and a fourth field which will depend on the type of entry. For QM atoms, the 
+# field will contain the element of the QM atom. For point charge lines, the field
+# will contain the charge of the point charge.
+# 
+# OUTPUT: The expected output file whould be placed in the same directory as the 
+# input file, and should be named "*inputfile*.result" (meaning it will have the
+# same path and name of the input file, plus the suffix '.result'). This file
+# should have, on its first line, the energy of the system, and on the following 
+# X lines (where X is the number of QM atoms passed in the input file), four (4) 
+# fields: the x, y and z components of the TOTAL FORCE applied on that atom, and 
+# on the fourth field, the charge of the atom. If the user indicates that charges
+# from the QM software should not be used (see "QMChargeMode"), the fourth field
+# should have zeroes, but should not be empty.
+
+
+#######  ORCA Files Format
+
+# ORCA input file is composed of a variable-sized header, where all info pertaining
+# to the simulation is specified, and a corrdinates section, where all atom positions
+# are specified, along with some informations regarding the overall system (namely
+# the charge and multiplicity of the system).
+#
+# We will get the header and the initial part of the coordinates section from 
+# lines specified in this script. The data would be gathered and writen by NAMD,
+# but in this case we need to previously prepare it from hardcoded values or by 
+# reading from another file.
+# 
+# ORCA expects a separate file with point charge information. This info is provided
+# by NAMD in the same data file, so we will extract it and paste it in a separate 
+# file in the same folder. The name and path to the point charge file needs to be 
+# included in the configuration portion of ORCA's input file.
+
+##################  Imports
+
+from sys import argv as sargv
+from sys import exit 
+from os.path import dirname
+import subprocess as sp
+
+################## Hardcoded parameters
+
+# Here we set up a template for the header and options of the input file.
+
+# The choice of method and basis-set, as well as parallelization and charge
+# output are all made here.
+# For ORCA, the keyword "ENGRAD" is essential for the proper output of gradients
+# over QM atoms.
+
+orcaConfigLines1 = """\
+!  B3LYP 6-31G* Grid4 PAL4 EnGrad TightSCF
+%output 
+  PrintLevel Mini 
+  Print [ P_Mulliken ] 1
+  Print [ P_AtCharges_M ] 1
+end
+
+"""
+
+# We set up the template for the point charge file indicator
+
+orcaConfigLines2 = "%pointcharges \""
+
+# And set up the header for the coordinates section of the input file.
+
+orcaConfigLines3 = """\
+%coords
+  CTyp xyz
+  Charge 0.000000
+  Mult 1.000000
+  Units Angs
+  coords
+
+"""
+
+################## Processing of file names
+
+inputFilename = sargv[1]
+
+#print(inputFilename)
+
+directory = dirname(inputFilename)
+
+# Prepares the name of the configuration file based on full path to data file provided by NAMD
+orcaInFileName = directory + "/"
+orcaInFileName += "qmmm.input"
+
+# Prepares the name of the point charge file based on full path to data file provided by NAMD
+#%pointcharges "/dev/shm/NAMD_test/0/qmmm_0.input.pntchrg"
+pcFileName = orcaInFileName
+pcFileName += ".pntchrg"
+
+orcaConfigLines2 += pcFileName + "\"\n"
+
+
+# Prepares the name of the output file based on full path to data file provided by NAMD.
+# This is where we will direct all ORCA output, so that we can grab the atom charge
+# information to pass back to NAMD.
+orcaOutFileName = orcaInFileName
+orcaOutFileName += ".TmpOut"
+
+# Prepares the name of the gradient file based on full path to data file provided by NAMD.
+# This is where ORCA will write the gradient information on QM atoms.
+orcaGradFileName = orcaInFileName
+orcaGradFileName += ".engrad"
+
+# Prepares the file name for the file which will be read by NAMD
+finalResFileName = inputFilename
+finalResFileName += ".result"
+
+#print("orcaInFileName:",orcaInFileName)
+#print("pcFileName:",pcFileName)
+#print("orcaOutFileName:",orcaOutFileName)
+#print("orcaGradFileName:",orcaGradFileName)
+#print("finalResFileName:",finalResFileName)
+
+################## Reading and parsing NAMD's data ; Writing ORCA's Input
+
+# Reads NAMD data
+infile = open(inputFilename,"r") 
+
+line = infile.readline()
+
+# Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
+numQMatms = int(line.split()[0])
+# Gets number of point charges
+numPntChr = int(line.split()[1].replace("\n",""))
+
+#print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
+
+# stores all lines written to ORCA's input file
+outLinesQM = []
+
+# stores all lines written to ORCA's point charge file
+outLinesPC = []
+
+# The first line in the point charge file is composed only of the total number
+# of point charges in the file.
+outLinesPC.append(str(numPntChr) + "\n")
+
+# Identation
+ident = "  "
+
+lineIndx = 1
+for line in infile:
+    
+    posx = line.split()[0]
+    posy = line.split()[1]
+    posz = line.split()[2]
+    
+    if lineIndx <= numQMatms:
+        
+        # ORCA's format requires the fileds to be ordered begining with the
+        # atom's element symbol, and followed by the XYZ coordinates.
+                
+        element = line.split()[3].replace("\n","")
+        
+        outLinesQM.append(ident + " ".join([element,posx,posy,posz]) + "\n")
+        
+    else:
+        
+        # ORCA's format requires the fileds to be ordered begining with the
+        # charge, and followed by the XYZ coordinates.
+        
+        charge = line.split()[3]
+        
+        outLinesPC.append(" ".join([charge,posx,posy,posz]) + "\n")
+    
+    lineIndx += 1
+    
+
+# Finalizes the formating for ORCA's input file, one "end" to terminate the
+# atomic coordinates, and nother to terminate the section.
+outLinesQM.append(ident + "end" + "\n")
+outLinesQM.append("end" + "\n")
+
+
+infile.close()
+
+###
+
+with open(orcaInFileName,"w") as outQMFile:
+    
+    outQMFile.write(orcaConfigLines1)
+    outQMFile.write(orcaConfigLines2)
+    outQMFile.write(orcaConfigLines3)
+    
+    for line in outLinesQM:
+        outQMFile.write(line)
+
+with open(pcFileName,"w") as outPCFile:
+    
+    for line in outLinesPC:
+        outPCFile.write(line)
+
+
+
+################## Run ORCA
+
+# We first move the shell to the target directory where calculations are to be
+# performed
+cmdline = "cd " + directory + "; "
+# Then we run orca with our output file receiving all standard output.
+cmdline += "/data/Programas/orca_3_0_3_linux_x86-64/orca "
+cmdline += orcaInFileName + " > " + orcaOutFileName
+
+#print("command:",cmdline)
+
+proc = sp.Popen(args=cmdline, shell=True)
+proc.wait()
+
+
+# Runs a secondary process
+cmdline2 = "/home/melomcr/Research/NAMD_QMMM/Scripts/saveDipole.sh "
+cmdline2 += orcaInFileName
+
+proc2 = sp.Popen(args=cmdline2, shell=True)
+
+################## Process ORCA results
+
+# Gets system energy and gradients from ORCA's output file.
+
+gradFile = open(orcaGradFileName,"r")
+    
+# Skips to the line with number of atoms
+for i in range(3):
+    gradFile.readline()
+
+orcaNumQMAtms = int(gradFile.readline().replace("\n",""))
+
+# Runs basic sanity check
+if orcaNumQMAtms != numQMatms:
+    print("ERROR: Expected",numQMatms,"but found",orcaNumQMAtms,"atoms in engrad file!")
+    exit(1)
+
+# Skips to the line with final system energy
+for i in range(3):
+    gradFile.readline()
+
+finalEnergy = gradFile.readline().replace("\n","").strip()
+
+print("ORCA energy: ", finalEnergy,"Eh")
+
+# All energies are given in Eh (Hartree)
+# NAMD needs energies in kcal/mol
+# The conversion factor is 627.509469
+finalEnergy = str( float(finalEnergy) * 627.509469 )
+
+print("ORCA energy: ", finalEnergy,"kcal/mol")
+
+# Skips to the lines with gradients
+for i in range(3):
+    gradFile.readline()
+
+# Stores all gradients (yes, it would be faster with numpy, but this is just an example)
+grads = []
+for i in range(orcaNumQMAtms):
+    
+    grads.append( list() )
+    
+    # All *GRADIENTS* are given in Eh/a0 (Hartree over Bohr radius)
+    # NAMD needs *FORCES* in kcal/mol/angstrons
+    # The conversion factor is -1*627.509469/0.529177 = -1185.82151
+    for j in range(3):
+        gradComp = gradFile.readline().replace("\n","").strip()
+        gradComp = float(gradComp) * -1185.82151
+        grads[i].append( str(gradComp) )
+
+gradFile.close()
+
+
+###
+
+# Gets atom charges from the temporary output file.
+tmpOutFile = open(orcaOutFileName,"r")
+
+qmCharges = []
+
+# Iterates ultil we find the section of output that contains atomic partial 
+# charges for QM atoms
+chargeSection = False
+
+iterate = True
+
+while iterate:
+    
+    line = tmpOutFile.readline()
+    
+    if line.find("MULLIKEN ATOMIC CHARGES") != -1:
+        chargeSection = True
+        
+        # Skips a dividing line
+        line = tmpOutFile.readline()
+        
+        continue
+    
+    if chargeSection:
+        qmCharges.append(line.split()[3].replace("\n","").strip())
+        pass
+    
+    if len(qmCharges) == numQMatms:
+        break
+
+tmpOutFile.close()
+
+# Writes the final output file that will be read by NAMD.
+
+finFile = open(finalResFileName,"w")
+
+# The first line constains the enegy
+finFile.write(finalEnergy + "\n")
+
+# And all follwoing lines contain gradients and partial charges
+for i in range(numQMatms):
+    
+    finFile.write(" ".join(grads[i]) + " " + qmCharges[i] + "\n")
+
+finFile.close()
+
+
+##########
+
+# Makes sure the secondary process has ended befire we return to NAMD.
+proc2.wait()
+
+
+exit(0)
diff --git a/lib/qmmm/run_gaussian.py b/lib/qmmm/run_gaussian.py
new file mode 100755 (executable)
index 0000000..df16275
--- /dev/null
@@ -0,0 +1,260 @@
+#!/usr/bin/python
+
+# by Maximilian Scheurer (mscheurer@ks.uiuc.edu), Dec 2016
+
+##################  Imports
+
+from sys import argv as sargv
+from sys import exit
+from os.path import dirname
+import subprocess as sp
+import re
+import numpy as np
+
+
+################## Command file parameters
+
+gaussianConfigLines1 = """\
+%NProcShare=4
+# B3LYP/6-31G* Charge Force
+
+qmmm
+
+0 1
+"""
+
+# if you specified an excited state ( TD(NStates=2,Root=1) ) in the com file and want to give its gradient to NAMD,
+# specify the number of the excited state here,
+# so that the python script can read the corresponding excited state energy
+# if no TD calculation is requested, set excitedState to 0
+excitedState=0
+
+
+################## Processing of file names
+gaussianWhitespace = "\n\n"
+
+inputFilename = sargv[1]
+
+#print(inputFilename)
+
+directory = dirname(inputFilename)
+
+# Prepares the name of the configuration file based on full path to data file provided by NAMD
+gaussianInFileName = directory + "/"
+gaussianInFileName += "qmmm.com"
+
+# Name of the Gaussian log file
+gaussianOutFileName = gaussianInFileName +".log"
+
+# Prepares the file name for the file which will be read by NAMD
+finalResFileName = inputFilename
+finalResFileName += ".result"
+
+################## Reading and parsing NAMD's data ; Writing gaussian's Input
+
+# Reads NAMD data
+infile = open(inputFilename,"r")
+
+line = infile.readline()
+
+# Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
+numQMatms = int(line.split()[0])
+# Gets number of point charges
+numPntChr = int(line.split()[1].replace("\n",""))
+
+# print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
+
+# stores all lines written to gaussian's input file
+outLinesQM = []
+
+# Identation
+ident = "  "
+
+lineIndx = 1
+pointChargeList = []
+pointChargeDict = {}
+charges = []
+for line in infile:
+
+    posx = line.split()[0]
+    posy = line.split()[1]
+    posz = line.split()[2]
+
+    if lineIndx <= numQMatms:
+
+        # gaussian's format requires the fileds to be ordered begining with the
+        # atom's element symbol, and followed by the XYZ coordinates.
+
+        element = line.split()[3].replace("\n","")
+
+        outLinesQM.append(" ".join([element,posx,posy,posz]) + "\n")
+
+    else:
+        #output linebreak to separate atoms from charges
+        if lineIndx == numQMatms+1:
+            outLinesQM.append("\n")
+
+        # gaussian's format requires the fields to be ordered begining with the
+        # XYZ, and followed by the charge .
+        pos = " ".join(line.split()[0:3])
+        charge = line.split()[3]
+        charges.append(charge)
+        if pos in pointChargeDict:
+               print "double occurence: ", pos, pointChargeDict[pos] , " with new charge ", charge
+               pointChargeDict[pos] += float(charge)
+               print "new merged charge: ", pointChargeDict[pos]
+        else:
+               pointChargeDict[pos] = float(charge)
+
+
+    lineIndx += 1
+
+cnp = np.array(charges,dtype=float)
+print "Sum of the charges: ", np.sum(cnp)
+
+for k in pointChargeDict:
+       c = pointChargeDict[k]
+       p = k.split()
+       pointChargeList.append([p[0],p[1],p[2],str(c)])
+       outLinesQM.append(" ".join([p[0],p[1],p[2],'{0:.16f}'.format(c)]) + "\n")
+
+# print len(pointChargeList)
+infile.close()
+
+###
+
+with open(gaussianInFileName,"w") as outQMFile:
+
+    outQMFile.write(gaussianConfigLines1)
+
+    for line in outLinesQM:
+        outQMFile.write(line)
+
+    outQMFile.write(gaussianWhitespace)
+
+################## Run gaussian
+
+# We first move the shell to the target directory where calculations are to be
+# performed
+cmdline = "cd " + directory + "; "
+# Then we run gaussian with our output file receiving all standard output.
+
+# we probably need to set some environment variables:
+import subprocess, os
+current_env = os.environ.copy()
+current_env["GAUSS_EXEDIR"] = "/path/to/gaussian/exedir"
+# subprocess.Popen(my_command, env=my_env)
+
+cmdline += "/path/to/gaussian/exedir/g09 "
+cmdline += gaussianInFileName + " " + gaussianOutFileName
+
+# print "command:", cmdline
+proc = sp.Popen(args=cmdline, shell=True, env=current_env)
+proc.wait()
+
+########### READING Gaussian output
+excitedStateEnergy=0
+
+tmpOutFile = open(gaussianOutFileName,"r")
+
+qmCharges = []
+gradient = []
+
+# Bohr radius for conversion
+a0 = 0.52917721067
+conversionHB = 627.509469/a0
+selfEnergyGaussian = 0.0
+
+# Iterates ultil we find the section of output that contains atomic partial
+# charges for QM atoms
+chargeSection = False
+
+gradientSection = False
+
+scfenergyFound = False
+
+excitedFound = False
+
+iterate = True
+
+selfEnergyFound = False
+
+while iterate:
+
+    line = tmpOutFile.readline()
+
+    if line.find("Mulliken charges:") != -1:
+        chargeSection = True
+        # Skips a dividing line
+        line = tmpOutFile.readline()
+        continue
+
+    if not scfenergyFound:
+        if line.find("SCF Done") != -1:
+            reg = re.search("SCF Done:(.+)=(.*)([-+]?\d*\.\d+|\d+)(.+)A.U.",line)
+            scfenergy = 627.509469 * float(reg.group(2).strip())
+            print "SCF energy: ", scfenergy
+            scfenergyFound = True
+
+    if not excitedFound and excitedState != 0:
+        if line.find("Excited State") != -1:
+            line = line.strip().replace(":","").split()
+            if int(line[2]) != excitedState:
+                continue
+            line =tmpOutFile.readline()
+            # check if we really requested the excited state that the gradient will be calculated for!
+            while line.find("This state for optimization and/or second-order correction.") == -1:
+                line = tmpOutFile.readline()
+            line = tmpOutFile.readline()
+            line = line.strip()
+            reg = re.search("([-+]?\d*\.\d+|\d+)(.+)",line)
+            excitedStateEnergy = 627.509469 * float(reg.group(1))
+            print "Excited State energy: " , excitedStateEnergy
+            excitedFound = True
+
+    if line.find("Center     Atomic                   Forces (Hartrees/Bohr)") != -1:
+        gradientSection = True
+        #skip 2 lines
+        line = tmpOutFile.readline()
+        line = tmpOutFile.readline()
+        line = tmpOutFile.readline()
+
+    if gradientSection:
+        line = line.strip().split()
+        # don't switch the sign of the number, as Gaussian prints F = -dE/dx, so we can directly pass the numbers to NAMD
+        gradient.append( [ str( float(line[2])*conversionHB ), str( float(line[3])*conversionHB ), str( float(line[4])*conversionHB ) ])
+        if len(gradient) == numQMatms:
+            gradientSection = False
+            iterate = False
+
+    if chargeSection:
+        qmCharges.append(line.split()[2].replace("\n","").strip())
+
+    if len(qmCharges) == numQMatms:
+        chargeSection = False
+
+    if selfEnergyFound != True and line.find("Self energy of the charges") != -1:
+       line = line.strip().split()
+       selfEnergyGaussian = float(line[-2])
+       print "Gaussian self energy of the charges: ", selfEnergyGaussian
+        selfEnergyFound = True
+
+tmpOutFile.close()
+
+finFile = open(finalResFileName,"w")
+
+if excitedState == 0:
+       finFile.write(str( scfenergy - selfEnergyGaussian*627.509469 ) + "\n")
+else:
+       finFile.write(str( excitedStateEnergy - selfEnergyGaussian*627.509469 ) + "\n")
+
+for i in range(numQMatms):
+
+    finFile.write(" ".join(gradient[i]) + " " + qmCharges[i] + "\n")
+
+finFile.close()
+
+
+##########
+
+exit(0)
diff --git a/lib/qmmm/run_qchem.py b/lib/qmmm/run_qchem.py
new file mode 100755 (executable)
index 0000000..d43bea3
--- /dev/null
@@ -0,0 +1,301 @@
+#!/usr/bin/python
+
+# by Maximilian Scheurer (mscheurer@ks.uiuc.edu), June 2017
+
+##################  Imports
+
+from sys import argv as sargv
+from sys import exit
+from os.path import dirname
+import subprocess as sp
+import re
+import numpy as np
+import math
+
+
+################## Command file parameters
+qchemNumberOfCores = 2
+
+qchemConfigLine1 = """\
+$rem
+JOBTYPE force
+METHOD b3lyp
+BASIS 6-31G*
+SCF_CONVERGENCE 7
+$end
+
+$molecule
+0 1
+"""
+
+
+
+
+
+# if you specified an excited state ( CIS, ADC or whatever ) in the input file and want to give its gradient to NAMD,
+# specify the number of the excited state here,
+# so that the python script can read the corresponding excited state energy
+# if no exc. state calculation is requested, set excitedState to 0
+excitedState=0
+
+# add the following lines to qchemConfigLine1 to calculate e.g. 4 excited states with TDDFT,
+# where the gradient will be calculated on the 1st excited state
+# CIS_N_ROOTS 4
+# CIS_STATE_DERIV 1
+# CIS_SINGLETS true
+# CIS_TRIPLETS false
+
+
+
+################## Processing of file names
+
+inputFilename = sargv[1]
+
+#print(inputFilename)
+
+directory = dirname(inputFilename)
+
+# Prepares the name of the configuration file based on full path to data file provided by NAMD
+qchemInFileName = directory + "/"
+qchemInFileName += "qmmm.in"
+
+# Name of the qchem log file
+qchemOutFileName = qchemInFileName +".out"
+
+# Prepares the file name for the file which will be read by NAMD
+finalResFileName = inputFilename
+finalResFileName += ".result"
+
+################## Reading and parsing NAMD's data ; Writing gaussian's Input
+
+# Reads NAMD data
+infile = open(inputFilename,"r")
+
+line = infile.readline()
+
+# Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
+numQMatms = int(line.split()[0])
+# Gets number of point charges
+numPntChr = int(line.split()[1].replace("\n",""))
+
+# print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
+
+# stores all lines written to gaussian's input file
+outLinesQM = []
+
+# Identation
+ident = "  "
+
+lineIndx = 1
+pointChargeList = []
+pointChargeDict = {}
+charges = []
+for line in infile:
+
+    posx = line.split()[0]
+    posy = line.split()[1]
+    posz = line.split()[2]
+
+    if lineIndx <= numQMatms:
+
+        # qchem's format requires the fileds to be ordered begining with the
+        # atom's element symbol, and followed by the XYZ coordinates.
+
+        element = line.split()[3].replace("\n","")
+
+        outLinesQM.append(" ".join([element,posx,posy,posz]) + "\n")
+
+    else:
+        #output linebreak to separate atoms from charges
+        if lineIndx == numQMatms+1:
+            outLinesQM.append("\n")
+
+        # qchems's format requires the fields to be ordered begining with the
+        # XYZ, and followed by the charge .
+        pos = " ".join(line.split()[0:3])
+        charge = line.split()[3]
+        charges.append(charge)
+        if pos in pointChargeDict:
+               print "double occurence: ", pos, pointChargeDict[pos] , " with new charge ", charge
+               pointChargeDict[pos] += float(charge)
+               print "new merged charge: ", pointChargeDict[pos]
+        else:
+               pointChargeDict[pos] = float(charge)
+
+
+    lineIndx += 1
+
+cnp = np.array(charges,dtype=float)
+print "Sum of the charges: ", np.sum(cnp)
+
+outLinesQM.append("$end \n \n $external_charges \n")
+for k in pointChargeDict:
+       c = pointChargeDict[k]
+       p = k.split()
+       pointChargeList.append([p[0],p[1],p[2],str(c)])
+       outLinesQM.append(" ".join([p[0],p[1],p[2],'{0:.16f}'.format(c)]) + "\n")
+
+outLinesQM.append("$end")
+# print len(pointChargeList)
+infile.close()
+
+###
+
+with open(qchemInFileName,"w") as outQMFile:
+
+    outQMFile.write(qchemConfigLine1)
+
+    for line in outLinesQM:
+        outQMFile.write(line)
+
+    # outQMFile.write(gaussianWhitespace)
+
+################## Run qchem
+
+# We first move the shell to the target directory where calculations are to be
+# performed
+cmdline = "cd " + directory + "; "
+# Then we run qchem with our output file receiving all standard output.
+
+# we probably need to set some environment variables:
+import subprocess, os
+current_env = os.environ.copy()
+current_env["QCSCRATCH"] = directory
+# subprocess.Popen(my_command, env=my_env)
+
+#load qchem from modules
+cmdline += "qchem -nt %d " % qchemNumberOfCores
+cmdline += qchemInFileName + " " + qchemOutFileName
+
+print "command:", cmdline
+proc = sp.Popen(args=cmdline, shell=True, env=current_env)
+proc.wait()
+
+########### READING QCHEM output
+excitedStateEnergy=0
+
+tmpOutFile = open(qchemOutFileName,"r")
+
+qmCharges = []
+gradient = []
+
+# Bohr radius for conversion
+a0 = 0.52917721067
+conversionHB = 627.509469/a0
+selfEnergy = 0.0
+
+# Iterates ultil we find the section of output that contains atomic partial
+# charges for QM atoms
+chargeSection = False
+
+gradientSection = False
+
+scfenergyFound = False
+
+excitedFound = False
+
+iterate = True
+
+selfEnergyFound = False
+
+while iterate:
+
+    line = tmpOutFile.readline()
+
+    if line.find("Ground-State Mulliken Net Atomic Charges") != -1:
+        chargeSection = True
+        # Skips 3 dividing lines
+        for x in range(3):
+            line = tmpOutFile.readline()
+        continue
+
+    if not scfenergyFound:
+        if line.find("SCF   energy in the final basis set") != -1:
+            scfenergy = 627.509469 * float(line.strip().split()[-1])
+            print "SCF energy: ", scfenergy
+            scfenergyFound = True
+
+    if not excitedFound and excitedState != 0:
+        if line.find("Excited State") != -1:
+            line = line.strip().replace(":","").split()
+            if int(line[2]) != excitedState:
+                continue
+            line =tmpOutFile.readline()
+            # check if we really requested the excited state that the gradient will be calculated for!
+            while line.find("This state for optimization and/or second-order correction.") == -1:
+                line = tmpOutFile.readline()
+            line = tmpOutFile.readline()
+            line = line.strip()
+            reg = re.search("([-+]?\d*\.\d+|\d+)(.+)",line)
+            excitedStateEnergy = 627.509469 * float(reg.group(1))
+            print "Excited State energy: " , excitedStateEnergy
+            excitedFound = True
+
+    if line.find("Gradient of") != -1:
+        gradientSection = True
+        # line = tmpOutFile.readline()
+
+
+    if gradientSection:
+        gradient = np.zeros(shape=(numQMatms,3))
+
+        ln = 0
+        blockColumnNumber = -1
+        atomLineIndex = 0
+        gradientRead = 0
+        while (gradientRead == 0):
+            line = tmpOutFile.readline()
+            line = line.strip().split()
+            if ln % 4 == 0:
+                blockColumnNumber = len(line)
+                # print(line, blockColumnNumber)
+                atomLineIndex += blockColumnNumber
+            else:
+                xyzIndex = (ln % 4)-1
+                # print("current coordinate: "+ str(xyzIndex))
+                currentColumnNumber = blockColumnNumber
+                columnCounter = 1
+                for atomIndex in range(atomLineIndex-currentColumnNumber, atomLineIndex):
+                    gradient[atomIndex,xyzIndex] = float(line[columnCounter])
+                    #print(atomIndex,columnCounter)
+                    if atomIndex == (numQMatms-1) and xyzIndex == 2:
+                        gradientRead = 1
+                    columnCounter += 1
+            ln+=1
+
+        iterate = 0
+
+    if chargeSection:
+        qmCharges.append(line.split()[2].replace("\n","").strip())
+
+    if len(qmCharges) == numQMatms:
+        chargeSection = False
+
+    if selfEnergyFound != True and line.find("Charge-charge energy") != -1:
+       line = line.strip().split()
+       selfEnergy = float(line[-2])
+       print "Self energy of the charges: ", selfEnergy
+        selfEnergyFound = True
+
+tmpOutFile.close()
+
+finFile = open(finalResFileName,"w")
+
+if excitedState == 0:
+       finFile.write(str( scfenergy - selfEnergy*627.509469 ) + "\n")
+       print "Corrected SCF energy: ", scfenergy-selfEnergy*627.509469
+else:
+       finFile.write(str( excitedStateEnergy - selfEnergy*627.509469 ) + "\n")
+
+forces=np.multiply(-1.0*conversionHB,gradient)
+forces=forces.tolist()
+#print forces
+for i in range(numQMatms):
+    finFile.write(' '.join('{0:.7f}'.format(c) for c in forces[i]) + " " + qmCharges[i] + "\n")
+
+finFile.close()
+
+
+##########
+
+exit(0)
diff --git a/lib/qmmm/run_terachem.py b/lib/qmmm/run_terachem.py
new file mode 100755 (executable)
index 0000000..734855f
--- /dev/null
@@ -0,0 +1,295 @@
+#!/usr/bin/python
+
+# by Maximilian Scheurer (mscheurer@ks.uiuc.edu), June 2017
+
+##################  Imports
+
+from sys import argv as sargv
+from sys import exit
+from os.path import dirname
+import subprocess as sp
+import numpy as np
+from scipy import spatial
+
+# TeraChem does not print the point-charge self energy, so we need to
+# calculate it ourselves. This may take a lot of time if many point charges
+# are included due to O(N^2) scaling!
+def computePointChargeSelfEnergy(pclist):
+      energy = 0.0
+      length = len(pclist)
+      pc = np.array(pclist,dtype=np.float64)
+      a = np.arange(0,length,1)
+      coords = pc[:,:3]
+      r = spatial.distance.cdist(coords,coords)/0.52917721067
+      idc = np.triu_indices(length,1)
+      energy = np.sum(pc[:,3][idc[0]]*pc[:,3][idc[1]]/(r[idc]))
+      return energy
+
+
+
+################## Hardcoded parameters
+
+# Here we set up a template for the header and options of the input file.
+
+excitedState=0
+
+tcConfigLines1 = """\
+basis 6-31G*
+coordinates qmmm.xyz
+pointcharges point_charges
+charge 0
+spinmult 1
+poptype mulliken
+method b3lyp
+dftd no
+run gradient
+end
+"""
+
+inputFilename = sargv[1]
+
+directory = dirname(inputFilename)
+
+# Prepares the name of the configuration file based on full path to data file provided by NAMD
+tcInFileName = directory + "/"
+tcInFileName += "qmmm_tc.in"
+
+coordFile = directory + "/qmmm.xyz"
+
+# Prepares the name of the point charge file based on full path to data file provided by NAMD
+#%pointcharges "/dev/shm/NAMD_test/0/qmmm_0.input.pntchrg"
+pcFileName = directory + "/point_charges"
+
+
+# Prepares the name of the output file based on full path to data file provided by NAMD.
+# This is where we will direct all tc output, so that we can grab the atom charge
+# information to pass back to NAMD.
+tcOutFileName = tcInFileName
+tcOutFileName += ".out"
+
+# Prepares the file name for the file which will be read by NAMD
+finalResFileName = inputFilename
+finalResFileName += ".result"
+
+# print("tcInFileName:",tcInFileName)
+# print("pcFileName:",pcFileName)
+# print("tcOutFileName:",tcOutFileName)
+#print("tcGradFileName:",tcGradFileName)
+# print("finalResFileName:",finalResFileName)
+
+################## Reading and parsing NAMD's data ; Writing tc's Input
+# Reads NAMD data
+infile = open(inputFilename,"r")
+
+line = infile.readline()
+
+# Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
+numQMatms = int(line.split()[0])
+# Gets number of point charges
+numPntChr = int(line.split()[1].replace("\n",""))
+
+#print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
+
+# stores all lines written to tc's input file
+outLinesQM = []
+
+# stores all lines written to tc's point charge file
+outLinesPC = []
+
+# The first line in the point charge file is composed only of the total number
+# of point charges in the file.
+
+# Identation
+ident = "  "
+
+lineIndx = 1
+lineIndx = 1
+pointChargeList = []
+pointChargeDict = {}
+charges = []
+for line in infile:
+
+    posx = line.split()[0]
+    posy = line.split()[1]
+    posz = line.split()[2]
+
+    if lineIndx <= numQMatms:
+
+
+        element = line.split()[3].replace("\n","").capitalize()
+
+        outLinesQM.append(ident + " ".join([element,posx,posy,posz]) + "\n")
+
+    else:
+
+
+        charge = line.split()[3]
+        pos = " ".join(line.split()[0:3])
+        charges.append(charge)
+        if pos in pointChargeDict:
+                print("double occurence: ", pos, pointChargeDict[pos] , " with new charge ", charge)
+                pointChargeDict[pos] += float(charge)
+                print("new merged charge: ", pointChargeDict[pos])
+        else:
+                pointChargeDict[pos] = float(charge)
+
+
+    lineIndx += 1
+
+
+cnp = np.array(charges,dtype=float)
+print("Sum of the charges: ", np.sum(cnp))
+
+outLinesPC.append(str(len(pointChargeDict)) + "\n\n")
+for k in pointChargeDict:
+        c = pointChargeDict[k]
+        p = k.split()
+        pointChargeList.append([p[0],p[1],p[2],str(c)])
+        outLinesPC.append(" ".join(['{0:.4f}'.format(c),p[0],p[1],p[2]]) + "\n")
+
+pcSelfEnergy = computePointChargeSelfEnergy(pointChargeList)
+print("self-energy: " , str(pcSelfEnergy) , " a.u.")
+
+infile.close()
+
+###
+
+with open(tcInFileName,"w") as outQMFile:
+
+    outQMFile.write(tcConfigLines1)
+
+
+with open(pcFileName,"w") as outPCFile:
+
+    for line in outLinesPC:
+        outPCFile.write(line)
+
+with open(coordFile,"w") as cFile:
+       cFile.write(str(numQMatms))
+       cFile.write("\n\n")
+       for line in outLinesQM:
+               cFile.write(line)
+
+################## Run tc
+
+# We first move the shell to the target directory where calculations are to be
+# performed
+import os
+current_env = os.environ.copy()
+
+cmdline = "cd " + directory + "; "
+# Then we run tc with our output file receiving all standard output.
+
+# here, scipy is not regularly in the python path :(
+cmdline += "$TeraChem/bin/terachem "
+cmdline += tcInFileName + " > " + tcOutFileName
+
+#print("command:",cmdline)
+
+proc = sp.Popen(args=cmdline, shell=True,env=current_env)
+proc.wait()
+
+
+print("terachem finished running.")
+################## Process tc results
+
+scfEnergy = 0
+iterate = True
+
+gradientSection = False
+excSection = False
+grads = []
+a0 = 0.52917721067
+conversionHB = 627.509469/a0
+eEnergy = 0
+resultFile = open("scr/results.dat","r")
+mainFile = open(tcOutFileName,"r")
+
+if not excitedState:
+       while iterate:
+               line = mainFile.readline()
+
+               if line.find("dE/dX") != -1:
+                       gradientSection = True
+                       line = mainFile.readline()
+
+               if gradientSection == True:
+                       g = line.strip().split()[0:3]
+                       grads.append([-1.0*float(x)*conversionHB for x in g])
+                       if len(grads) == numQMatms:
+                               gradientSection = False
+                               break
+
+
+iterate = True
+while iterate:
+       line = resultFile.readline()
+       if line == None:
+               break
+       if line.find("Ground state energy (a.u.):") != -1:
+               line = resultFile.readline()
+               scfEnergy = (float(line.strip()) - pcSelfEnergy)*627.509469
+               if excitedState == 0:
+                       break
+
+       if line.find("dE/dX") != -1:
+               gradientSection = True
+               line = resultFile.readline()
+
+       if gradientSection == True:
+               g = line.strip().split()[1:4]
+               grads.append([-1.0*float(x)*conversionHB for x in g])
+               if len(grads) == numQMatms:
+                       gradientSection = False
+                       if excitedState == 0:
+                               break
+
+       if line.find("Root      Energy") != -1:
+               line = resultFile.readline()
+               excSection = True
+
+       if excSection == True:
+               state, energy = line.strip().split()
+               if int(state) == excitedState:
+                       # eV to kcal/mol
+                       eEnergy = float(energy) * 23.06054819
+                       break
+
+
+resultFile.close()
+finalEnergy = scfEnergy + eEnergy
+print(scfEnergy,eEnergy,finalEnergy)
+###
+
+# Gets atom charges from the charge_mull.xls output file.
+mullikenFile = open("scr/charge_mull.xls","r")
+
+qmCharges = []
+
+iterate = True
+
+while iterate:
+    line = mullikenFile.readline()
+    qmCharges.append(line.split()[2].replace("\n","").strip())
+
+    if len(qmCharges) == numQMatms:
+        break
+
+mullikenFile.close()
+# Writes the final output file that will be read by NAMD.
+finFile = open(finalResFileName,"w")
+
+# The first line constains the enegy
+finFile.write(str(finalEnergy) + "\n")
+
+# And all follwoing lines contain gradients and partial charges
+for i in range(numQMatms):
+
+    finFile.write(" ".join(str(x) for x in grads[i]) + " " + qmCharges[i] + "\n")
+
+finFile.close()
+
+
+##########
+
+exit(0)
diff --git a/lib/qmmm/saveDipole.sh b/lib/qmmm/saveDipole.sh
new file mode 100755 (executable)
index 0000000..433724e
--- /dev/null
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+targDir=`dirname "$1"`
+ID=`basename "$targDir"`
+
+echo "Input file: $1"
+echo "Target Dir: $targDir"
+echo "Target ID: $ID"
+
+dp=`cat $1.TmpOut | grep "Total Dipole Moment" | sed -e 's/ \{1,\}/,/g' -e 's/^.*:,\(.*\)/\1/'`
+echo "Dipole moment: $dp"
+
+if [ ! -f "/home/melomcr/Research/NAMD_QMMM/TestSystem/dipoleMoments_$ID.csv" ] ; then
+echo "x,y,z" > /home/melomcr/Research/NAMD_QMMM/TestSystem/dipoleMoments_$ID.csv
+fi
+
+echo "$dp" >> /home/melomcr/Research/NAMD_QMMM/TestSystem/dipoleMoments_$ID.csv
\ No newline at end of file
index 1140c55..8222da7 100644 (file)
@@ -30,7 +30,7 @@ ug_dynamics.tex ug_forcefield.tex ug_performance.tex ug_userdef.tex \
 ug_intro.tex ug_io.tex ug_macros.tex ug_runit.tex ug_sample.tex ug_start.tex \
 ug_xplor.tex ug_fenergy.tex ug_psfgen.tex psfgen_macros.tex ug_alchemy.tex \
 ug_colvars.bib ug_colvars_macros.tex \
-ug_colvars.tex ug_analysis.tex ug_gbis.tex ug_constantph.tex
+ug_colvars.tex ug_analysis.tex ug_gbis.tex ug_constantph.tex ug_qmmm.tex
 
 PDFFIGS = figures/fmaOn.pdf figures/pairlistdist.pdf \
 figures/shifting.pdf figures/switching.pdf figures/dual_top.pdf
diff --git a/ug/figures/ChargeGroupDiagram.png b/ug/figures/ChargeGroupDiagram.png
new file mode 100644 (file)
index 0000000..32c475f
Binary files /dev/null and b/ug/figures/ChargeGroupDiagram.png differ
diff --git a/ug/figures/MultipleQMDiagram.png b/ug/figures/MultipleQMDiagram.png
new file mode 100644 (file)
index 0000000..6c9f211
Binary files /dev/null and b/ug/figures/MultipleQMDiagram.png differ
diff --git a/ug/figures/OptionsDiagram.png b/ug/figures/OptionsDiagram.png
new file mode 100644 (file)
index 0000000..715caef
Binary files /dev/null and b/ug/figures/OptionsDiagram.png differ
diff --git a/ug/figures/QMMMBond_2bonds-01.png b/ug/figures/QMMMBond_2bonds-01.png
new file mode 100644 (file)
index 0000000..3efacde
Binary files /dev/null and b/ug/figures/QMMMBond_2bonds-01.png differ
diff --git a/ug/figures/hybrid_qmmm_diagram.png b/ug/figures/hybrid_qmmm_diagram.png
new file mode 100644 (file)
index 0000000..c5eca8f
Binary files /dev/null and b/ug/figures/hybrid_qmmm_diagram.png differ
index 6028f63..92d64d1 100644 (file)
--- a/ug/ug.bib
+++ b/ug/ug.bib
 @string{JPCA = "J.~Phys.\ Chem.~A"}
 @string{JPCB = "J.~Phys.\ Chem.~B"}
 @string{JPCL = "J.~Phys.\ Chem.\ Lett."}
-@string{JCTC = "J.~Chem.\ Theor.\ Comp."}
+@string{JCTC = "J.~Chem.\ Theory\ Comput."}
 @string{CPL = "Chem.\ Phys.\ Lett."}
 @string{CP = "Chem.\ Phys."}
 
+@string{ACIEE = "Angew.\ Chem.\ Int.\ Ed.\ Engl."}
+@string{NATMET = "Nat.\ Methods"}
+@string{SR = "Sci.\ Rep."}
+@string{JCAMD = "J.~Comp.-Aided\ Mol.\ Design"}
+@string{WIRCMS = "Wiley Interdiscip.\ Rev.:\ Comput.\ Mol.\ Sci."}
+
 @inproceedings{CRUZ93,
   author={Carolina Cruz-Neira and Dan J. Sandin and Thomas A. DeFanti},
   title={Unknown},
@@ -1181,3 +1187,79 @@ URL = {http://www.eurekaselect.com/node/75602/article}
 
 % end REST2 references
 
+% begin QM/MM references
+
+@article{SENN2009,
+  title={QM/MM methods for biomolecular systems},
+  author={Senn, Hans Martin and Thiel, Walter},
+  journal=ACIEE,
+  volume={48},
+  number={7},
+  pages={1198--1229},
+  year={2009},
+  publisher={Wiley Online Library}
+}
+
+@article{MELO2018,
+   author = {Marcelo Melo and Rafael Bernardi and Till Rudack and Maximilian Scheurer and Christoph Riplinger and James Phillips and Julio Maia and Gerd Rocha and J. Ribeiro and John Stone and Frank Nesse and Klaus Schulten and Zaida Luthey-Schulten},
+   title = {{NAMD} goes quantum: An integrative suite for {QM/MM} simulations},
+   journal = NATMET,
+   volume={15},
+   pages={351--354},
+   year={2018},
+   doi = {},
+   tbreference={720},
+   tbstatus={Published.},
+   note={},
+   pmcid={PMC6095686},
+   annote={}
+}
+
+@article{RIBE2016,
+author={Joao V. Ribeiro and Rafael C. Bernardi and Till Rudack and John E. Stone and James C. Phillips and Peter L. Freddolino and Klaus Schulten},
+title={Qwik{MD}-Integrative Molecular Dynamics Toolkit for Novices and Experts},
+journal=SR,
+volume={6},
+pages={26536},
+year={2016},
+note={},
+doi={10.1038/srep26536},
+tbreference={687},
+tbstatus={Published.},
+pmcid={PMC4877583},
+annote={}
+}
+
+@article{STEW90,
+  title={MOPAC: a semiempirical molecular orbital program},
+  author={Stewart, James JP},
+  journal=JCAMD,
+  volume={4},
+  number={1},
+  pages={1--103},
+  year={1990},
+  publisher={Springer}
+}
+
+@article{MAIA2012,
+  title={GPU linear algebra libraries and GPGPU programming for accelerating MOPAC semiempirical quantum chemistry calculations},
+  author={Maia, Julio Daniel Carvalho and Urquiza Carvalho, Gabriel Aires and Mangueira Jr, Carlos Peixoto and Santana, Sidney Ramos and Cabral, Lucidio Anjos Formiga and Rocha, Gerd B},
+  journal=JCTC,
+  volume={8},
+  number={9},
+  pages={3072--3081},
+  year={2012},
+  publisher={ACS Publications}
+}
+
+@article{NEES2012,
+  author={Frank Neese},
+  journal=WIRCMS,
+  volume={2},
+  year=2012,
+  pages={73--78},
+  title={The {ORCA} program system}
+}
+
+% end QM/MM references
+
index 39befb4..1f7977f 100644 (file)
--- a/ug/ug.tex
+++ b/ug/ug.tex
 \newpage
 \input{ug_constantph}
 
+% QM/MM
+\newpage
+\input{ug_qmmm}
+
 % Runtime analysis
 \newpage
 \input{ug_analysis}
diff --git a/ug/ug_qmmm.tex b/ug/ug_qmmm.tex
new file mode 100644 (file)
index 0000000..4c36e3d
--- /dev/null
@@ -0,0 +1,841 @@
+% QM/MM
+% Documentation copied from https://www.ks.uiuc.edu/Research/qmmm/.
+
+\section{Hybrid QM/MM Simulations}
+\label{section:qmmm}
+
+% This brief introduction aims at providing some basic context 
+% to the following description of capabilities and commands 
+% available in NAMD's QM/MM interface.
+
+Even though molecular mechanics (MM) force-fields are based on
+quantum mechanical calculations and experimental observations,
+only quantum mechanics (QM) can give a complete and accurate
+understanding of many biochemical processes, particularly those
+involving chemical reactions or charge redistribution.
+Nevertheless, even with the advanced hardware technology available today,
+the computational cost of studying nanosecond-long dynamics of entire
+systems relying solely on QM methodologies is usually prohibitive.
+A common route to circumvent this cost barrier is to confine the
+QM formalism to a sub-region of a system and to include the effects
+of the surrounding system through MM simulations,
+leading to hybrid QM/MM simulations~\cite{SENN2009}.
+
+NAMD's comprehensive QM/MM suite~\cite{MELO2018} was developed to provide
+easy setup, visualization and analysis of QM/MM simulations through
+the graphical user interface VMD/QwikMD~\cite{RIBE2016},
+and a broad range of QM methods through NAMD's new ``QMForces" module.
+The QM/MM interface in NAMD supports the simulation of many
+independent QM regions, and smooth integration with a vast collection
+of enhanced sampling methods. In hybrid QM/MM simulations,
+NAMD offloads part of its standard force and energy calculations
+to a QM program, either through native interfaces to
+MOPAC~\cite{STEW90,MAIA2012} or ORCA~\cite{NEES2012},
+or through a flexible generic interface requiring a wrapper script,
+where exemplary Python wrappers are provided for Gaussian,
+TeraChem and Q-CHEM.
+Multiple QM-MM coupling schemes are implemented, allowing for both
+mechanically and electrostatically embedded QM regions to be used
+(see description in Nature Methods~\cite{MELO2018}).
+QM/MM simulations require the same input files used for classical MD,
+with additional options in the configuration file.
+QM and MM atoms covalently bound are usually treated by redistributing
+the MM atom's charge over its nearest MM neighbors and by capping
+the QM atom with a hydrogen atom,
+as shown in Figure~\ref{fig:hybrid_qmmm} for a solvated
+tri-alanine QM/MM calculation using the NAMD/ORCA interface.
+Tests of the QM/MM interface for accuracy, stability and performance,
+are provided as supporting information in Nature Methods~\cite{MELO2018}. 
+
+If employing NAMD QM/MM please cite:
+\begin{quote}
+NAMD goes quantum: An integrative suite for hybrid simulations.
+Melo*, M. C. R.; Bernardi*, R. C.; Rudack T.; Scheurer, M.; Riplinger, C.;
+Phillips, J. C.; Maia, J. D. C.; Rocha, G. D.; Ribeiro, J. V.;
+Stone, J. E.; Neese, F.; Schulten, K.; Luthey-Schulten, Z.;
+Nature Methods, 2018 (doi:10.1038/nmeth.4638)
+\end{quote}
+
+\begin{figure}[tbp]
+\centering
+\includegraphics[width=6in]{figures/hybrid_qmmm_diagram.png}
+\caption[Hybrid QM/MM NAMD]{%
+Graphical representation of NAMD-ORCA interconnection.
+Only the contribution of MM charges beyond rmax are
+calculated by NAMD (via PME), with the direct electrostatic
+calculation performed by ORCA.
+The image assumes the charge shift redistribution scheme,
+where the partial charge of the linking MM atom is shifted
+to its nearest MM neighbors.
+}
+\label{fig:hybrid_qmmm}
+\end{figure}
+
+
+\subsection{Division of Labor}
+
+The basic idea behind a hybrid QM/MM simulation in NAMD is to use 
+a classical force field to treat the classical atoms in the system 
+(or ``MM atoms"), and pass the information that describes the quantum 
+atoms in the system (or ``QM atoms") to a Quantum Chemistry (QC) software, 
+which is expected to produce gradients for all QM atoms, as well as 
+the total energy of the QM region (and optionally partial charges). 
+All bonded and non-bonded interactions among MM atoms are handled 
+by NAMD's force field. Similarly, all bonded and non-bonded interactions 
+among QM atoms are handled by the QC software in its chosen theory level. 
+Treatment of covalent bonds between QM and MM atoms will be described 
+in a following section.
+
+The non-bonded interactions between QM and MM atoms are handled differently, 
+and can be modified and regulated by the user. 
+Van der Waals interactions are always calculated, and can be done 
+using either the default force field parameters, or specific 
+(user-defined) parameters for QM atoms. 
+Parameter modifications for QM atoms have been proposed 
+in order to compensate for over-polarization that these atoms 
+may exhibit in hybrid QM/MM simulations. 
+Larger van der Waals radii and/or shallower well depths 
+should then be provided for all element types that occur 
+among QM atoms (see the ``qmVdwParams" keyword).
+
+
+\subsection{Mechanical and Electrostatic Embedding}
+
+Electrostatic interactions between QM and MM atoms deserve 
+a more detailed discussion due to the abundance and diversity 
+of available alternatives. 
+The first decision to be made is whether there will be electrostatic 
+interactions between the two portions of a system, QM and MM. 
+In the ``mechanical embedding" scheme, only positions and elements 
+of atoms in the QM region are passed on to the chosen QC software 
+for energy and force calculations. 
+This way, QM and MM atoms share only van der  Waals interactions.
+
+% Figure 1 for QM/MM
+\begin{figure}[tbp]
+\centering
+\includegraphics[width=5in]{figures/OptionsDiagram.png}
+\caption[Diagram of classical point charge options.]{%
+Diagram of options that control the use and manipulation of 
+classical point charges. Default values are indicated below their 
+respective keyword. ``Cutoff" and ``SwitchDist" are keywords used in 
+NAMD to configure the calculations of electrostatic and van der Waals 
+interactions. 
+}
+\label{fig:qmmm_options}
+\end{figure}
+
+In the ``electrostatic embedding" scheme, on the other hand, 
+the partial charges of MM atoms surrounding all QM atoms 
+are used to approximate the electrostatic environment where QM atoms 
+are found (the scheme is selected with the ``qmElecEmbed" keyword). 
+See Figure~\ref{fig:qmmm_options}.
+This process can be customized in a variety 
+of ways, the first of which is deciding if a smoothing function will 
+be used to avoid an abrupt decay in electrostatic force 
+due to the cutoff used in the selection of surrounding point charges
+(this option is activated with the ``qmSwitching" keyword).
+
+Classical point charge utilization can be further customized 
+by choosing which smoothing function will be used, 
+and if the total charge of selected partial charges should be 
+modified to (A) have a whole charge or (B) have a complementary 
+charge to that of the QM region, so that the sum of charges 
+from QM atoms and classical partial charges add to zero
+(see Figure~\ref{fig:qmmm_options}).
+
+With electrostatic embedding, QM atoms are influenced by the charges 
+in the classical region. In order to balance the forces acting on the 
+system, NAMD uses partial charges for the QM atoms to calculate 
+the electrostatic interaction with classical point charges. 
+There are two possibilities for the origin of the QM partial charges: 
+the original partial charges found in the force field parameter files 
+can be used, or updated partial charges can be gathered at each step 
+from the QC software output 
+(controllable through the ``qmChargeMode" keyword). 
+The continuous update in charge distribution allows for a partial 
+re-parameterization of the targeted molecule at each time step, 
+which can lead to an improved description of the interactions 
+of a ligand as it repositions over the surface of a protein, 
+for example, or as it moves through a membrane.
+
+In case PME is activated by the user, NAMD will automatically apply 
+the necessary corrections to handle the QM region, allowing it to be 
+influenced by long range interactions from the entire system.
+
+
+\subsection{Covalent Bonds Divided by the QM/MM Barrier}
+
+Hybrid QM/MM simulations of biomolecular systems often present 
+situations where only a portion of a molecule should be treated 
+quantum mechanically, usually to save computational resources 
+since the cost of simulating QM regions rises rapidly with the 
+number of simulated toms. In order to deal with chemical bonds 
+that are split by the QM/MM division of the biomolecular system, 
+that is, bonds that have one atom in the quantum (QM) region and 
+another in the classical (MM) region (we will call these ``QM/MM bonds"), 
+NAMD makes approximations to the molecular system in order to bridge 
+differences in simulation type (QM vs.\ MM), and minimize errors involved 
+in the QM/MM division of the system
+(Figure~\ref{fig:qmmm_bond_treatment} A and B).
+
+\begin{figure}[tbp]
+\centering
+\includegraphics[width=5in]{figures/QMMMBond_2bonds-01.png}
+\caption[Treatment of QM/MM bonds]{%
+Treatment of QM/MM bonds.
+A) Illustration of all atoms in the vicinity of the QM/MM bond, 
+colored by element: cyan for carbon, white for hydrogen, 
+blue for nitrogen and red for oxygen. 
+B) To the left, in blue, is the region that will be treated with 
+the chosen QC software. To the right, in red, the region treated 
+classically by NAMD. The bond in gray is the one crossing the 
+QM/MM division. The atom marked as \textbf{QM1} is the quantum atom 
+directly connected to the classical atom on the other side of the 
+QM/MM division. Analogously, the atom marked as \textbf{MM1} is the 
+classical atom directly connected to the quantum atom on the other 
+side of the QM/MM division. Atoms marked as \textbf{MM2} are directly 
+bonded to the \textbf{MM1} atom, and atoms marked \textbf{MM3} are 
+directly bonded to \textbf{MM2} atoms. 
+C) \textbf{Z1} method. Ignored partial charges are indicated in the image 
+with a gray sphere taking the place of its respective classical atom. 
+Directly adjacent to MM1 is a green sphere representing the link atom 
+that is placed along the QM1-MM1 covalent bond. All remaining partial 
+charges representing classical atoms that are passed on to the QC software
+are indicated in purple spheres. 
+D) \textbf{Z2} method.  
+E) \textbf{Z3} method. 
+F) RCD method. Virtual point charges, are represented in yellow spheres. 
+The text indicates the total charge placed at each position, 
+where \textbf{q} indicates the charge of the MM1 atom and \textbf{q2} 
+represents the partial charge of the MM2 atom at that position. 
+The yellow spheres at MM2 atom positions indicate their partial charge 
+has been changed from its original value. 
+G) CS method. Since in this case the virtual point charges are placed 
+very close to the MM2 atom position, the yellow spheres representing 
+them show significant overlapping.
+}
+\label{fig:qmmm_bond_treatment}
+\end{figure}
+
+
+\subsubsection{Link Atoms}
+
+As previously mentioned, the information regarding atoms in the QM region 
+is passed on to the chosen QC software, that is, their respective 
+positions and element types, but in order to maintain (or approximate) 
+the effect of the chemical bond between the QM atom and the MM atom, 
+NAMD creates and places a link atom (usually a hydrogen) along the 
+``broken" QM/MM bond. The user can fine-tune this process by choosing 
+the method of placement of the link atom and even the element of 
+such atom (keywords ``qmBondDist" and ``qmLinkElement").
+
+The introduction of the link atom will invariably place it very near 
+the classical atom involved in the QM/MM bond, therefore the use 
+and placement of partial charges from classical atoms becomes
+highly relevant. Under the mechanical embedding scheme, the QC software 
+only receives the atoms in the QM region and the link atoms created
+to approximate QM/MM bonds, so no manipulation of partial charges is 
+required. On the other hand, usual QM/MM simulations are done under 
+the electrostatic embedding scheme, in which case the partial charges
+of classical atoms involved in the QM/MM bonds and classical atoms 
+directly connected to them require special treatment.
+
+
+\subsubsection{Point Charge Alterations}
+
+Several methods have been proposed to handle this situation, and the 
+QM/MM interface developed here implements the most widely accepted ones.
+One can be chosen using the ``qmBondScheme" keyword
+(Figure~\ref{fig:qmmm_bond_treatment} C to G). 
+In all implemented methods, the classical atom participating in the 
+QM/MM bond (MM1 atom) does not have its partial charge passed on 
+to the QC software, since this would create excessive repulsion 
+(or attraction) on the link atom. This is, in fact, the entirety of 
+the ``Z1" method: ignoring the partial charge of the MM1 atom. 
+Analogously, ``Z2" and ``Z3" ignore all partial charges up to 
+MM2 and MM3 atoms, respectively
+(Figure~\ref{fig:qmmm_bond_treatment} C to E).
+
+The Redistributed Charge and Dipole (RCD) method
+(Figure~\ref{fig:qmmm_bond_treatment} F)
+is more elaborate, as it rearranges the partial charge of the 
+MM1 atom (indicated as $q$) so that the total charge of the region 
+is maintained as well as the dipole moments of the bonds between 
+MM1 and MM2 atoms. This is done by creating ``virtual" point charges, 
+which are passed on to the QC software as if they represented 
+partial charges of classical atoms. More specifically, the RCD method 
+creates a virtual point charge in the middle of all MM1-MM2 bonds 
+with a charge of $2q/n$, where $n$ is the number of 
+MM2 atoms connected to MM1, and also subtracts a charge $q/n$ 
+from each of the MM2 atoms, so that the total charge of the region 
+remains constant while approximating the dipole moment of the 
+MM1-MM2 bonds. This way there will be no point charge placed at the 
+position of the MM1 atom, but its partial charge is not simply removed, 
+it is redistributed.
+
+A similar approach is taken by the Charge Shifting (CS) method 
+(Figure~\ref{fig:qmmm_bond_treatment} G).
+In this case, the MM1 partial charge is equally 
+distributed across the MM2 atoms, and two virtual point charges 
+are placed along the direction of the MM1-MM2 bond, one before 
+the MM2 atom and one after, each one with a charge of $+q/n$ 
+and $-q/n$ respectively. This method will also keep the 
+total charge of the region constant while trying to preserve 
+the local dipoles formed by all MM1-MM2 bonds.
+
+
+\subsubsection{Link Atom Charge and Charge Groups}
+
+Along with the gradient over all QM atoms, NAMD can also use 
+the partial charge derived from the QC calculation to update 
+the charge distribution of atoms in the QM region. 
+When a QM/MM bond exists, however, part of the charge of the 
+region will be placed on the link atom, and in order to keep the 
+charge of the QM region constant, the link atom charge is 
+re-distributed on the QM region. This seemingly simple mechanism 
+can cause problems unless special care is be taken when deciding 
+which bond will mark the division of QM and MM regions.
+
+Many force fields divide the topologies of biomolecules in 
+``charge groups"
+(Figure~\ref{fig:qmmm_chargegroups} A and B).
+What this means is that 
+not only will the partial charges of all atoms of a molecule 
+add up to the whole number that represents the charge of the molecule, 
+they will also add up to whole numbers in sub groups of atoms 
+(look for the ``GROUP" statements in
+\url{http://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/node24.html}
+to see an example). 
+Therefore, one needs to make sure that the chosen QM/MM bond(s) sits
+in between ``charge groups", so the total sum of partial charges 
+of atoms defining a QM region is a whole number. This is especially 
+important in order to keep the total charge of the system constant. 
+Since the QC calculation will always distribute a whole charge over 
+all atoms of a system (QM atoms plus a link atom), if the partial charge 
+of QM atoms is not initially a whole number, it will be forced into 
+a whole number after the first QC step, where the charge of the link atom 
+is distributed over the QM region. This will create a mismatch between 
+QM and MM charges, changing the total charge of the entire system 
+(QM plus MM regions).
+
+\begin{figure}[tbp]
+\centering
+\includegraphics[width=5in]{figures/ChargeGroupDiagram.png}
+\caption[Charge Groups and QM/MM Bonds]{%
+Charge Groups and QM/MM Bonds.
+A) Illustration of aspartate and the distribution of charge over its 
+atoms as in CHARMM36 force field parameters. Circles in red indicate 
+oxygen atoms, blue indicate nitrogen atoms, cyan for carbon atoms, 
+and white for hydrogen atoms. ``Bond 1" indicates the peptide bond,
+``Bond 2" indicates the one between the alpha carbon and the peptide bond 
+nitrogen, and ``Bond 3" the bond between the alpha carbon and the peptide 
+bond carbon. 
+B) Charge groups are indicated with dashed squares, 
+along with their total charges. 
+C) Depiction of the atoms in the QM region if Bonds 1 and 3 are used to 
+separate it from the MM region. The total charge of QM region is $-1$. 
+D) Depiction of QM region if the same is defined by Bonds 2 and 3. 
+In this case, the total charge of QM region is $-1.16$.
+}
+\label{fig:qmmm_chargegroups}
+\end{figure}
+
+An example can be seen in Figure~\ref{fig:qmmm_chargegroups},
+bonds 1 and 3 are chosen as the QM/MM bonds,
+the charge distribution seen in Figure~\ref{fig:qmmm_chargegroups} C shows 
+a whole charge for the QM region (and consequently for the MM region). 
+Therefore, any charge placed on link atoms can be redistributed 
+to the QM atoms with no change in total system charge. However, 
+if bonds 2 and 3 are chosen for the QM/MM bond
+(Figure~\ref{fig:qmmm_chargegroups} D), 
+the charge of the MM region would be $+1.16$, while the charge 
+of the QM region would be $-1.16$. Since the QC calculation would 
+place a pre-determined whole charge on the region ($-1$, in this case), 
+the updated charge of the QM region after the first simulation step 
+would change the total charge of the system to $+0.16$, in this example.
+
+
+\subsection{Custom Quantum Chemistry Software}
+
+In order to offer the broad range of tools and technologies present 
+in NAMD to all researchers who develop and/or employ specialized 
+Quantum Chemistry tools, the QM/MM interface is prepared to utilize 
+any QC tool that can be wrapped in a script that converts input and 
+output files to specified formats. This flexible interface will 
+improve development and testing of new tools, as well as their 
+quick integration utilization in hybrid dynamics.
+
+We currently provide 
+in the \texttt{libs/qmmm/} directory
+(and mirrored at
+\url{http://www.ks.uiuc.edu/Research/qmmm/Scripts/})
+Python wrapper scripts for GAUSSIAN, TeraChem, and Q-Chem.
+Other wrapper scripts can be generated, based on these templates,
+in any other language,
+as long as they are provided to NAMD in an executable form.
+Although natively supported,
+we also provide a python wrapper script for ORCA,
+with extended comments explaining the format in which NAMD will write 
+data for the QC software and the format in which NAMD expects to find 
+the results.
+
+
+\subsection{Independent QM Regions}
+
+Aiming at large macromolecular simulations that could take advantage 
+of localized QM resolution, NAMD allows the user to set up multiple 
+independent QM regions in the same molecular system. For example, 
+one could study a multimeric complex that contains several active sites 
+and have all active sites be calculated with a chosen QC software 
+simultaneously (Figure~\ref{fig:qmmm_multiple_grid}).
+Each active site would be calculated 
+independently of all others, by its own execution of the QC software, 
+keeping the calculation cost low and without impacting the overall 
+efficiency of the simulation, since all QM regions would be calculated 
+in parallel.
+
+\begin{figure}[tbp]
+\centering
+\includegraphics[width=5in]{figures/MultipleQMDiagram.png}
+\caption[Diagram of Multiple Grid Regions]{%
+Diagram of Multiple QM Regions. 
+The illustration depicts a hetero-hexameric complex 
+(light and dark green triangles) that combine to create three active sites 
+(in gray). Active sites are bound to a target molecule (purple triangle). 
+The inner and outer dashed circles represent, respectively, the boundary 
+of a QM region and the limits of the classical point charge shell around 
+that region.
+}
+\label{fig:qmmm_multiple_grid}
+\end{figure}
+
+Identifying the different QM regions and which atoms belong to each 
+one of them can be simply accomplished in the input PDB file, 
+or in a dedicated PDB file (keyword ``qmParamPDB"). Since each region 
+can contain different sets of molecules, their charges and multiplicities 
+are indicated separately (see keywords ``qmCharge" and ``qmMult").
+
+For simulations of large systems that are distributed across several 
+computer nodes, one can control how many independent QM regions are 
+calculated in each node. This would prevent large simulations from 
+running out of memory if two or more large QM regions are placed 
+in the same node (see keyword ``qmSimsPerNode").
+
+
+\subsection{Keywords}
+
+\begin{itemize}
+\setlength{\itemsep}{0.4cm}
+
+\item
+\NAMDCONFWDEF{qmForces}{Calculate QM?}{yes or no}{no}{%
+Turns on or off the QM calculations.
+}
+
+\item
+\NAMDCONF{qmParamPDB}{Set QM atoms}{PDB file}{%
+Name of an optional secondary PDB file where the OCCupancy
+or BETA column has the indications for QM or MM atoms.
+QM atoms should have an integer bigger than zero (0) and 
+MM atoms should have zero as the beta or occupancy field. 
+The same file may have indications for bonds between a QM atom 
+and an MM atom. This should be provided in case the PDB file 
+passed to the ``coordinates" keyword already has data on its 
+beta or occupancy columns, such as when a SMD simulations 
+is being performed.
+}
+
+\item
+\NAMDCONF{qmColumn}{Which column?}{``beta" or ``occ"}{%
+Indicates which column has the QM/MM field.
+Required.
+}
+
+\item
+\NAMDCONFWDEF{qmSimsPerNode}{Sims per node}{postive integer}{1}{%
+Number of independent simultaneous QM simulations per node.
+}
+
+\item
+\NAMDCONF{qmBondColumn}{Which bond column?}{``beta" or ``occ"}{%
+Indicates which column has the QM/MM bond information.
+This will tell NAMD which atoms are at the ends of a covalent bond
+split by the QM/MM barrier, with one atom being quantum and
+one being classical.
+There is no default value.
+If this parameter is provided,
+NAMD will parse options for dealing with QM/MM bonds.
+}
+
+\item
+\NAMDCONFWDEF{qmBondDist}{Use qmBondColumn value for distance?}%
+{``on'' or ``off''}{off}{%
+Indicates whether the value in the BondColumn will be used to define 
+the distance between the QM atom and the Link Atom that will replace 
+the MM atom in the QM system.
+}
+
+\item
+\NAMDCONFWDEF{qmBondValueType}{Does qmBondColumn value give length or ratio?}%
+{``len'' or ``ratio''}{len}{%
+Indicates if the values in the BondColumn represent either 
+the length (``len'') between the QM and link atoms or 
+the ratio (``ratio'') between the QM-MM distance and 
+the one which will be used as the QM-Link distance.
+}
+
+\item
+\NAMDCONFWDEF{qmLinkElement}{Set link atom element}%
+{string, for example ``4 9 Cl"}{H}{%
+User defined link atom element. Two syntaxes are allowed:
+if there is only one QM-MM bond, a string with the element symbol 
+is allowed (such as ``H" or ``Cl"). If there are two or more bonds, 
+the string needs to have the two atoms that compose the bond, 
+and then the element (such as ``4 9 Cl"). 
+The default element for all link atoms is hydrogen.
+}
+
+\item
+\NAMDCONFWDEF{qmBondScheme}{Select QM-MM bond scheme}%
+{``CS" or ``RCD" or ``Z1" or ``Z2" or ``Z3"}{CS}{%
+Indicates what will be the treatment given to QM-MM bonds in terms 
+of charge distribution and link atom creation and placement. 
+CS: Charge Shift Scheme; RCD: Redistributed Charge and Dipole method; 
+Z1: Only ignored MM1 partial charge, with no charge redistribution; 
+Z2: Ignores MM1 and all MM2 partial charges, with no charge redistribution; 
+Z3: Ignores MM1 and all MM2 and MM3 partial charges, 
+with no charge redistribution.
+}
+
+\item
+\NAMDCONFWDEF{qmElecEmbed}{Should point charges be used in QM?}%
+{``on" or ``off"}{on}{%
+Indicates if classical point charges should be used in QM calculations.
+}
+
+\item
+\NAMDCONFWDEF{qmSwitching}{Use switching on point charges?}%
+{``on" or ``off"}{off}{%
+This will scale down the point charges representing the classical 
+system as to replicate the switching procedure that NAMD applies 
+to all charged interaction (see ``switching").
+}
+
+\item
+\NAMDCONFWDEF{qmSwitchingType}{Set functional form of switching}%
+{``switch" or ``shift"}{shift}{%
+This option is used to decide which kind of function will be used 
+to scale down point charges sent to QM calculations. 
+SHIFT: This will "shift down" the entire shell of point charges 
+so that electrostactic interactions reach zero at the cutoff distance. 
+SWITCH: This will only change point charges in the sub-volume between 
+the switchdist and cutoff distance, so that electrostactic interactions 
+reach zero at the cutoff distance.
+}
+
+\item
+\NAMDCONFWDEF{qmPointChargeScheme}{Set point charge scheme}%
+{``none" or ``round" or ``zero"}{none}{%
+This option allows the user to decide if and how the point charges 
+presented to the QM system will be altered. NONE: Nothing will be done. 
+ROUND: This will change the most distant point charges so that the 
+total sum of point charges is a whole number. ZERO: This will adjust 
+the most distant point charges so that the total sum of point charges 
+is 0.
+}
+
+\item
+\NAMDCONF{qmBaseDir}{Set directory for file I/O}{directory path}{%
+This should be a fast read/write location, such as a RAM drive
+(/dev/shm on most linux distros). The user needs to make sure this 
+directory exists in the node(s) running the QM calculation(s).
+}
+
+\item
+\NAMDCONFWDEF{qmReplaceAll}{Use only QM gradients for forces?}%
+{``on" or ``off"}{off}{%
+Indicates to NAMD that ALL forces form NAMD will be ignored and 
+only the gradients from the QM software will be applied on the atoms. 
+This IS NOT NECESSARY in any regular QM/MM simulation, and will prevent 
+the use of any other feature from NAMD such as SMD.
+}
+
+\item
+\NAMDCONFWDEF{qmVdwParams}{Modify type names for QM atoms?}%
+{``on" or ``off"}{off}{%
+The QM code will change all QM atoms' van der Waals types to "q"+element 
+(e.g., all carbons will be qC and all hydrogens will be qH)
+for VdW interactions. This means that a parameter file with 
+epsilon and sigma values for these  atom types must be provided 
+along with the regular parameter files. For example, if using 
+CHARMM force field, the new file should be in CHARMM format.
+}
+
+\item
+\NAMDCONFWDEF{qmPCStride}{Set stride for point charge}{integer}{1}{%
+Sets a stride for new point charge determination. The same set of 
+classical atoms will be sent to QM calculations as point charges, 
+but with updated positions.
+}
+
+\item
+\NAMDCONFWDEF{qmCustomPCSelection}{Provide custom point charge selection?}%
+{``on" or ``off"}{off}{%
+Indicates that one or more file(s) will be provided with a custom 
+selection of point charges. Each file will have a selection for a single
+QM group. This selection will be kept during the entire simulation.
+}
+
+\item
+\NAMDCONF{qmCustomPCFile}{File for custom point charge selection}{PDB file}{%
+The file will have, in the ``qmColumn", the same QM ID provided for
+a single QM group. All other groups will have zero (0) in this column.
+In the second column (beta or occupancy), the classical or quantum atoms
+(from other QM regions) that need to be passed as point charges will be
+identified by a non-zero number.
+
+\vspace{-0.25em}
+\textbf{Example/Format:}
+\texttt{%
+\vspace{-1em}
+\small{
+\begin{tabbing}
+qmCustomPCFile \= system/system.customPC.1.pdb \\
+qmCustomPCFile \> system/system.customPC.2.pdb \\
+qmCustomPCFile \> system/system.customPC.3.pdb \\
+qmCustomPCFile \> system/system.customPC.4.pdb
+\end{tabbing}
+} % \small
+} % \texttt
+} % \NAMDCONF
+
+\item
+\NAMDCONFWDEF{qmLiveSolventSel}{Keep track of solvent?}{``on" or ``off"}{off}{%
+With Live Solvent Selection (LSS), NAMD will automatically keep track
+of the solvent molecules for all QM Groups, and will exchange classical
+solvent molecules with QM solvent molecules every "QMLSSFreq" steps.
+}
+
+\item
+\NAMDCONFWDEF{qmLSSResname}{Set residue name for LSS}{residue name}{TIP3}{%
+Indicates which residue name will be used in LSS.
+}
+
+\item
+\NAMDCONFWDEF{qmLSSFreq}{Set frequency of LSS}%
+{integer, multiple of stepspercycle}{100}{%
+Frequency of LSS. Must be a multiple of stepspercycle.
+}
+
+\item
+\NAMDCONFWDEF{qmLSSMode}{How solvent molecules are selected}%
+{``dist" or ``COM"}{dist}{%
+For LSS, this indicates how solvent molecules are selected.
+In all cases, the closest solvent molecules are selected,
+and if a classical solvent molecule is closer than a QM one,
+they are swaped. DIST: This mode will use the smallest distance
+between a solvent atom and a non-solvent QM atom to sort solvent
+molecules. This is best used when the non-solvent QM atoms form
+irregular volumes (when the COM is not very representatve),
+and/or volumes with high solvent accessibility (such as a drug,
+or a small peptide, in solution). COM: This mode will sort solvent
+molecules based on Center Of Mass distance between the solvent COM
+and the COM of a selection for each QM group (see ``qmLSSRef" keyword).
+Best used with small QM regions that have limited solvent accessibility,
+such as an active site.
+}
+
+\item
+\NAMDCONF{qmLSSRef}{Which residues for COM of QM atoms?}{string}{%
+This will indicate which residues are to be used in the determination
+of the COM of non-solvent QM atoms. Only these atoms will be used to
+determine the closest set of solvent molecules. The keyword takes a
+string composed of the QM group ID, the segment name and the residue ID.
+
+\vspace{-0.25em}
+\textbf{Example/Format:}
+\texttt{%
+\vspace{-1em}
+\small{
+\begin{tabbing}
+qmLSSRef \= "1 RP1 9" \\
+qmLSSRef \> "2 RP1 3" \\
+qmLSSRef \> "2 RP1 2" \\
+qmLSSRef \> "3 AP1 9" \\
+qmLSSRef \> "3 AP1 3" \\
+qmLSSRef \> "4 AP1 9"
+\end{tabbing}
+} % \small
+} % \texttt
+} % \NAMDCONF
+
+\item
+\NAMDCONF{qmConfigLine}{Pass string to QM configuration}{string}{%
+The string passed to "qmConfigLine" will be copied and pasted at the
+very begining of the configuration file for the chosen QM software
+if either ORCA or MOPAC are selected.
+
+\vspace{-0.25em}
+\textbf{Example/Format (QM/MM NAMD-ORCA):}
+\texttt{%
+\vspace{-1em}
+\small{
+\begin{tabbing}
+qmConfigLine \= "! PM3 ENGRAD" \\
+qmConfigLine \> "\%\%output PrintLevel Mini Print\textbackslash{[}P\_Mulliken\textbackslash{]} 1 Print\textbackslash{[}P\_AtCharges\_M\textbackslash{]} 1 end"
+\end{tabbing}
+} % \small
+} % \texttt
+
+\vspace{-1em}
+\textbf{Example/Format (QM/MM NAMD-MOPAC):}
+\texttt{%
+\vspace{-1em}
+\small{
+\begin{tabbing}
+qmConfigLine \= "PM7 XYZ T=2M 1SCF MOZYME CUTOFF=9.0 AUX LET GRAD QMMM GEO-OK" \\
+qmConfigLine \> "Test System"
+\end{tabbing}
+} % \small
+} % \texttt
+
+} % \NAMDCONF
+
+\item
+\NAMDCONF{qmMult}{Set multiplicity of QM region}{string}{%
+Multiplicity of the QM region. This is needed for proper construction
+of ORCA's input file. Each string must be composed of the QM region ID
+and its multiplicity.
+
+\vspace{-0.25em}
+\textbf{Example/Format:}
+\texttt{%
+\vspace{-1em}
+\small{
+\begin{tabbing}
+qmMult \= "1 1" \\
+qmMult \> "2 1" \\
+qmMult \> "3 1" \\
+qmMult \> "4 1"
+\end{tabbing}
+} % \small
+} % \texttt
+} % \NAMDCONF
+
+\item
+\NAMDCONF{qmCharge}{Set charge of each QM region}{string}{%
+Indicates the charge of each QM region. If no charge is provided
+for a QM region, NAMD calculates the total charge automatically
+based on the given parameter set. Each string must be composed
+of the QM region ID and its total charge.
+
+\vspace{-0.25em}
+\textbf{Example/Format:}
+\texttt{%
+\vspace{-1em}
+\small{
+\begin{tabbing}
+qmCharge \= "1 1" \\
+qmCharge \> "2 -1" \\
+qmCharge \> "3 1" \\
+qmCharge \> "4 -1"
+\end{tabbing}
+} % \small
+} % \texttt
+} % \NAMDCONF
+
+\item
+\NAMDCONF{qmSoftware}{Which QM software?}%
+{``mopac" or ``orca" or ``custom"}{%
+Required for QM/MM, this
+indicates which QM software should be used. In case the user wants
+to apply another QM code, this can be done by using the "custom"
+qmSoftware. In this case, NAMD will call the executable defined in
+the qmExecPath variable and will give it one argument: the full path
+to the input file. 
+
+INPUT: This input file will contain on the first line the number of
+QM atoms (X) and the number of point charges in the file (Y, which
+may be 0 or more), separated  by a space. The following X+Y lines
+will have four (4) fields: X, Y and Z coordinates, and a fourth field
+which will depend on the type of entry. For QM atoms, the  field will
+contain the element of the QM atom. For point charge lines, the field
+will contain the charge of the point charge.
+
+OUTPUT: The expected output file whould be placed in the same directory
+as the input file, and should be named \texttt{"*inputfile*.result"}
+(meaning it will have the same path and name of the input file, plus the
+suffix \texttt{".result"}). This file should have, on its first line,
+the energy of the system, and on the following X lines (where X is the
+number of QM atoms passed in the input file), four (4) fields: the x, y,
+and z components of the TOTAL FORCE applied on that atom, and on the
+fourth field, the charge of the atom. If the user indicates that charges
+from the QM software should not be used (see ``qmChargeMode"), the fourth
+field should have zeroes, but should not be empty. Energy should be
+in Kcal/mol and forces in Kcal/mol/Angstrom.
+}
+
+\item
+\NAMDCONF{qmExecPath}{Set path to QM executable}{path}{%
+Required for QM/MM, this
+indicates the path to the QM code executable.
+}
+
+\item
+\NAMDCONF{qmSecProc}{Set path to secondary executable}{path}{%
+Indicates a secondary executable that NAMD will call AFTER each
+QM software execution for each QM group. The executable is called
+with two arguments: the complete path and name of the input file
+used for each QM software execution; and the simulation step. This
+option can be used for an extra-processing at every step, e.g.,
+for saving all QM outputs every step.
+}
+
+\item
+\NAMDCONF{qmPrepProc}{Set path to initial executable}{path}{%
+Indicates an executable that must be called BEFORE the FIRST QM
+software execution of each QM group. The executable is called with
+one argument: the complete path and name of the input file used
+for each QM software execution. This can be used to setup a charge
+distribution for a molecule of interest, for example.
+}
+
+\item
+\NAMDCONFWDEF{qmChargeMode}{Set charge calculation mode}%
+{``none" or ``mulliken" or ``chelpg"}{mulliken}{%
+Charge calculation mode expected from the QM software. This
+indicates if charges should be read from the QM software and
+updated at every step, or if the original force field atom
+charges should be used. In case you are using ORCA, two charge
+options are allowed, Mulliken or CHELPG. We must know the kind
+of charge requested by the user so that the proper format is expected,
+read and processed. NONE: No charges are read from the QM software
+output and the original force field charges are preserved.
+MULLIKEN: This is the only other option for MOPAC and one possibility
+for ORCA. In case you are using the custom QM software interface,
+choose this option in order to use the charges calculated in the
+custom QM software, irrespective of the actual theory used for
+that calculation. CHELPG: This is a second possibility for ORCA.
+}
+
+\item
+\NAMDCONFWDEF{qmOutStride}{Set frequency of QM charge output}{integer}%
+{0 (not saving)}{%
+Frequency of QM charge output. A dedicated DCD file will be created
+to store the charge data for all QM atoms in the system. This
+independent frequency allows the user to store whole-system data
+at a larger stride to save time and space.
+}
+
+\item
+\NAMDCONFWDEF{qmPositionOutStride}{Set frequency of QM-only position output}%
+{integer}{0 (not saving)}{%
+Frequency of QM-only position output. A dedicated DCD file will be
+created to store the position data for all QM atoms in the system.
+This independent frequency allows the user to store whole-system
+data at a larger stride to save time and space.
+}
+
+\end{itemize}
+