modified on 18 February 2011 at 14:43 ••• 5,820 views

Qmmm PES analysis

From NWChem

Jump to: navigation, search

Contents

QM/MM Optimization and Transition States

Potential Energy Surface Analysis: Optimization | Transition States | Hessians and Frequency | NEB |


QM/MM optimization is based on multi-region optimization methodology and is invoked by

task qmmm <qmtheory> optimize

The overall algorithm involves alternating optimizations of QM and MM regions until convergence is achieved. This type of approach offers substantial savings compared to direct optimization of the entire system as a whole. In the simplest case of two regions (QM and MM) the algorithm is comprised of the following steps:

  1. Optimization of the QM region keeping MM region fixed
  2. Calculation of reduced electrostatic representation for the QM region (e.g. ESP charges)
  3. Optimization of MM region keeping QM region fixed
  4. Repeat from Step 1 until converged

The optimization process is controlled by the following keywords:

  • region - required keyword which specifies which regions will optimized and in which order.
  • maxiter - number of optimizations steps for each region within single optimization pass
  • ncycles - number of optimization cycles
  • density - electrostatic representation of the QM region during MM optimization
  • xyz - output of xyz structure files
  • convergence - convergence criteria


Here is an example QM/MM block that provides practical illustration of all these keywords for a generic optimization case where QM molecule(s) are embedded in the solvent

qmmm
  region  qm   solvent
  maxiter 10   3000
  ncycles 5
  density espfit
  xyz     foo
end

We have two regions in the system "qm" and "solvent" and we would like to optimize them both, thus the line

 region  qm   solvent

Our QM region is presumably small and the maximum number of iterations (within a single optimization pass) is set to 10. The solvent region is typically much larger (thousands of atoms) and the maximum number of iterations is set to a much large number 3000:

maxiter 10   3000

We would like to perform a total of 5 optimization passes, giving us a total of 5*10=50 optimization steps for QM region and 5*3000=15000 optimization steps for solvent region:

 ncycles 5

We are requesting QM region to be represented by point ESP charges during the solvent optimization:

density espfit

Finally we are requesting that the coordinates of the first region to be saved in the form of numbered xyz files:

xyz     foo

Example of QM/MM optimization

The example below illustrates QM/MM optimization at DFT/B3LYP level of theory for quantum water molecule embedded into 20 angstrom box of classical SPCE/E water molecules.

The restart (wtr_ref.rst) and topology (wtr.top) files are assumed to be generated elsewhere.

  start wtr

  permanent_dir ./perm
  scratch_dir ./data

  md
  system wtr_ref
  end

  basis
  * library "6-31G"
  end

  dft
  xc b3lyp
  end

  qmmm
  region  qm   solvent
  maxiter 10   1000
  ncycles 5
  density espfit
  xyz    foo
  end

  task qmmm dft optimize

QM/MM Hessian and Frequency Calculations

Potential Energy Surface Analysis: Optimization | Transition States | Hessians and Frequency | NEB |


Setup

QM/MM hessian and frequency calculations are invoked though the following task directives

task qmmm <qmtheory> hessian

or

task qmmm <qmtheory> freq

Only numerical implementation are supported at this point and will be used even in the absence of "numerical" keyword. Other than standard QM/MM directives no additional QM/MM input is required for default hessian/frequency for all the QM atoms. Using region keyword(first region only), hessian/frequency calculations can also be performed for classical solute atoms. If only classical atoms are involved density keyword can be utilized to enable frozen density or ESP charge representation for fixed QM region. Hessian/frequency calculations for solvent are not possible.

Examples

Example of QM/MM frequency calculation for classical region

This example illustrates QM/MM frequency calculation for Ala-Ser-Ala system. In this case instead of default QM region (see prepare block), the calculation is performed on classical solute part of the system as defined by region directive in QM/MM block. The electrostatic field from fixed QM region is represented by point ESP charges (see density directive). These ESP charges are calculated from wavefunction generated as a result of energy calculation.

 memory total 800 Mb

 start asa

 permanent_dir ./perm
 scratch_dir ./data

 #this will generate topology file (asa.top), restart (asa_ref.rst), and pdb (asa_ref.pdb) files.
 prepare
   source asa.pdb
   new_top new_seq
   new_rst
   modify atom 2:_CB quantum
   modify atom 2:2HB quantum
   modify atom 2:3HB quantum
   modify atom 2:_OG quantum
   modify atom 2:_HG quantum
   center
   orient
   solvate
   update lists
   ignore
   write asa_ref.rst
   write asa_ref.pdb   # Write out PDB file to check structure
 end
 task prepare

 md
   system asa_ref
 end

 basis "ao basis"
   * library "6-31G*"
 end

 dft
  print low
  iterations 500
 end

 qmmm
 region mm_solute
 density espfit
 end

 # run energy calculation to generate wavefunction file for subsequent ESP charge generation
 task qmmm dft energy
 task qmmm dft freq