

(7 intermediate revisions not shown) 
Line 1: 
Line 1: 
  [[Category:ReleaseDoc60]]
 +  #REDIRECT [[Release66:MP2]] 
  =MP2=
 +  
   +  
  There are (at least) three algorithms within NWChem that compute the MøllerPlesset (or manybody) perturbation theory secondorder correction<ref>Møller, C. and Plesset, M.S. (1934) "Note on an approximation treatment for manyelectron systems", ''Physical Review'' '''46''' 618622, doi:[http://dx.doi.org/10.1103/PhysRev.46.618 10.1103/PhysRev.46.618].</ref> to the HartreeFock energy (MP2). They vary in capability, the size of system that can be treated and use of other approximations
 +  
   +  
  * Semidirect  this is recommended for most large applications (up to about 2800 basis functions), especially on the IBM SP and other machines with significant disk I/O capability. Partially transformed integrals are stored on disk, multipassing as necessary. RHF and UHF references may be treated including computation of analytic derivatives. This is selected by specifying mp2 on the task directive, e.g.
 +  
   +  
  TASK MP2
 +  
   +  
  * Fullydirect<ref>Wong, A.T.; Harrison, R.J. and Rendell, A.P. (1996) "Parallel direct fourindex transformations", ''Theoretica Chimica Acta'' '''93''' 317331, doi:[http://dx.doi.org/10.1007/BF01129213 10.1007/BF01129213]</ref>  this is of utility if only limited I/O resources are available (up to about 2800 functions). Only RHF references and energies are available. This is selected by specifying direct_mp2 on the task directive, e.g.
 +  
   +  
  TASK DIRECT_MP2
 +  
   +  
  * Resolution of the identity (RI) approximation MP2 (RIMP2)<ref>Bernholdt, D.E. and Harrison, R.J. (1996) "Largescale correlated electronic structure calculations: the RIMP2 method on parallel computers", ''Chemical Physics Letters'' '''250''' (56) 477484, doi:[http://dx.doi.org/10.1016/00092614(96)000541 10.1016/00092614(96)000541]</ref>  this uses the RI approximation and is therefore only exact in the limit of a complete fitting basis. However, with some care, high accuracy may be obtained with relatively modest fitting basis sets. An RIMP2 calculation can cost over 40 times less than the corresponding exact MP2 calculation. RHF and UHF references with only energies are available. This is selected by specifying rimp2 on the task directive, e.g.,
 +  
   +  
  TASK RIMP2
 +  
   +  
  All three MP2 tasks share the same input block.
 +  
   +  
  MP2
 +  
  [FREEZE [[core] (atomic  <integer nfzc default 0>)] \
 +  
  [virtual <integer nfzv default 0>]]
 +  
  [TIGHT]
 +  
  [PRINT]
 +  
  [NOPRINT]
 +  
  [VECTORS <string filename default scfoutputvectors> \
 +  
  [swap [(alphabeta)] <integer pairlist>] ]
 +  
  [RIAPPROX <string riapprox default V>]
 +  
  [FILE3C <string filename default $file_prefix$.mo3cint">]
 +  
  [SCRATCHDISK <integer>]
 +  
  END
 +  
   +  
   +  
  ==FREEZE  Freezing orbitals==
 +  
   +  
  All MP2 modules support frozen core orbitals, however, only the direct MP2 and RIMP2 modules support frozen virtual orbitals.
 +  
   +  
  By default, no orbitals are frozen. The atomic keyword causes orbitals to be frozen according to the rules in the table below. Note that no orbitals are frozen on atoms on which the nuclear charge has been modified either by the user or due to the presence of an ECP. The actual input would be
 +  
   +  
  freeze atomic
 +  
   +  
  For example, in a calculation on <math>Si(OH)_2</math>, by default the lowest seven orbitals would be frozen (the oxygen 1<math>s</math>, and the silicon 1<math>s</math>, 2<math>s</math> and 2<math>p</math>).
 +  
   +  
   +  
  <center>
 +  
  { border='1' cellspacing='0'
 +  
  +Number of orbitals considered "core" in the "freeze by atoms" algorithm
 +  
   Period
 +  
   Elements
 +  
   Core Orbitals
 +  
   Number of Core
 +  
  
 +  
   0
 +  
   H  He
 +  
   
 +  
   0
 +  
  
 +  
   1
 +  
   Li  Ne
 +  
   1<math>s</math>
 +  
   1
 +  
  
 +  
   2
 +  
   Na  Ar
 +  
   1<math>s</math>2<math>s</math>2<math>p</math>
 +  
   5
 +  
  
 +  
   3
 +  
   K  Kr
 +  
   1<math>s</math>2<math>s</math>2<math>p</math>3<math>s</math>3<math>p</math>
 +  
   9
 +  
  
 +  
   4
 +  
   Rb  Xe
 +  
   1<math>s</math>2<math>s</math>2<math>p</math>3<math>s</math>3<math>p</math>4<math>s</math>3<math>d</math>4<math>p</math>
 +  
   18
 +  
  
 +  
   5
 +  
   Cs  Rn
 +  
   1<math>s</math>2<math>s</math>2<math>p</math>3<math>s</math>3<math>p</math>4<math>s</math>3<math>d</math>4<math>p</math>5<math>s</math>4<math>d</math>5<math>p</math>
 +  
   27
 +  
  
 +  
   6
 +  
   Fr  Lr
 +  
   1<math>s</math>2<math>s</math>2<math>p</math>3<math>s</math>3<math>p</math>4<math>s</math>3<math>d</math>4<math>p</math>5<math>s</math>4<math>d</math>5<math>p</math>6<math>s</math>4<math>f</math>5<math>d</math>6<math>p</math>
 +  
   43
 +  
  
 +  
  }
 +  
  </center>
 +  
   +  
   +  
  Caution: The rule for freezing orbitals "by atoms" are rather unsophisticated: the number of orbitals to be frozen is computed from the Table 16.1 by summing the number of core orbitals in each atom present. The corresponding number of lowestenergy orbitals are frozen  if for some reason the actual core orbitals are not the lowest lying, then correct results will not be obtained. From limited experience, it seems that special attention should be paid to systems including third and higher period atoms.
 +  
   +  
  The user may also specify the number of orbitals to be frozen by atom. Following the <math>Si(OH)_2</math> example, the user could specify
 +  
   +  
  freeze atomic O 1 Si 3
 +  
   +  
  In this case only the lowest four orbitals would be frozen. If the user does not specify the orbitals by atom, the rules default to Table 16.1.
 +  
   +  
  Caution: The system does not check for a valid number of orbitals per atom. If the user specifies to freeze more orbitals then are available for the atom, the system will not catch the error. The user must specify a logical number of orbitals to be frozen for the atom.
 +  
   +  
  The FREEZE directive may also be used to specify the number of core orbitals to freeze. For instance, to freeze the first 10 orbitals
 +  
   +  
  freeze 10
 +  
   +  
  or equivalently, using the optional keyword core
 +  
   +  
  freeze core 10
 +  
   +  
  Again, note that if the 10 orbitals to be frozen do not correspond to the first 10 orbitals, then the swap keyword of the VECTORS directive must be used to order the input orbitals correctly ([[#VECTORS  MO vectorsMO vectors]]).
 +  
   +  
  To freeze the highest virtual orbitals, use the virtual keyword. For instance, to freeze the top 5 virtuals
 +  
   +  
  freeze virtual 5
 +  
   +  
  Again, note that this only works for the directMP2 and RIMP2 energy codes.
 +  
   +  
  ==TIGHT  Increased precision==
 +  
   +  
  The TIGHT directive can be used to increase the precision in the MP2 energy and gradients.
 +  
   +  
  By default the MP2 gradient package should compute energies accurate to better than a microHartree, and gradients accurate to about five decimal places (atomic units). However, if there is significant linear dependence in the basis set the precision might not be this good. Also, for computing very accurate geometries or numerical frequencies, greater precision may be desirable.
 +  
   +  
  This option increases the precision to which both the SCF (from <math>10^{6}</math> to <math>10^{8}</math>) and CPHF (from <math>10^{4}</math> to <math>10^{6}</math>) are solved, and also tightens thresholds for computation of the AO and MO integrals (from <math>10^{9}</math> to <math>10^{11}</math>) within the MP2 code.
 +  
   +  
  ==SCRATCHDISK  Limiting I/O usage==
 +  
   +  
  This directive  used only in the semidirect algorithm  allows to limit the per process disk usage. Mandatory argument for this keyword is the maximum number of MBytes. For example, the following input line
 +  
   +  
  scratchdisk 512
 +  
   +  
  puts an upper limit of 512 MBytes to the semidirect MP2 usage of disk (again, on a per process base).
 +  
   +  
  ==PRINT and NOPRINT==
 +  
   +  
  The standard print control options are recognized. The list of recognized names are given in the table below.
 +  
   +  
   +  
  <center>
 +  
  { border='1' cellspacing='0'
 +  
  +Printable items in the MP2 modules and their default print levels
 +  
   Item
 +  
   Print Level
 +  
   Description
 +  
  
 +  
  
 +  
  
 +  
  
 +  
  
 +  
   '''RIMP2'''
 +  
  
 +  
  
 +  
  
 +  
   
 +  
  
 +  
  
 +  
   "2/3 ints"
 +  
   debug
 +  
   Partial 3center integrals
 +  
  
 +  
   "3c ints"
 +  
   debug
 +  
   MO 3center integrals
 +  
  
 +  
   "4c ints b"
 +  
   debug
 +  
   "B" matrix with approx. 4c integrals
 +  
  
 +  
   "4c ints"
 +  
   debug
 +  
   Approximate 4center integrals
 +  
  
 +  
   "amplitudes"
 +  
   debug
 +  
   "B" matrix with denominators
 +  
  
 +  
   "basis"
 +  
   high
 +  
  
 +  
  
 +  
   "fit xf"
 +  
   debug
 +  
   Transformation for fitting basis
 +  
  
 +  
   "geombas"
 +  
   debug
 +  
   Detailed basis map info
 +  
  
 +  
   "geometry"
 +  
   high
 +  
  
 +  
  
 +  
   "information"
 +  
   low
 +  
   General information about calc.
 +  
  
 +  
   "integral i/o"
 +  
   high
 +  
   File size information
 +  
  
 +  
   "mo ints"
 +  
   debug
 +  
  
 +  
  
 +  
   "pair energies"
 +  
   debug
 +  
   (working only in direct_mp2)
 +  
  
 +  
   "partial pair energies"
 +  
   debug
 +  
   Pair energy matrix each time it is updated
 +  
  
 +  
   "progress reports"
 +  
   default
 +  
   Report completion of timeconsuming steps
 +  
  
 +  
   "reference"
 +  
   high
 +  
   Details about reference wavefunction
 +  
  
 +  
   "warnings"
 +  
   low
 +  
   Nonfatal warnings
 +  
  }
 +  
  </center>
 +  
   +  
   +  
  ==VECTORS  MO vectors==
 +  
   +  
  All of the (supported) MP2 modules require use of converged canonical SCF (RHF or UHF) orbitals for correct results. The vectors are by default obtained from the preceding SCF calculation, but it is possible to specify a different source using the VECTORS directive. For instance, to obtain vectors from the file /tmp/h2o.movecs, use the directive
 +  
   +  
  vectors /tmp/h2o.movecs
 +  
   +  
  As noted above ([[#FREEZE  Freezing orbitalsFREEZE]]) if the SCF orbitals are not in the correct order, it is necessary to permute the input orbitals using the swap keyword of the VECTORS directive. For instance, if it is desired to freeze a total six orbitals corresponding to the SCF orbitals 15, and 7, it is necessary to swap orbital 7 into the 6th position. This is accomplished by
 +  
   +  
  vectors swap 6 7
 +  
   +  
  The swap capability is examined in more detail in [[SCF#VECTORS  input/output of MO vectorsInput/output of MO vectors]].
 +  
   +  
  ==RIMP2 fitting basis==
 +  
   +  
  The RIMP2 method requires a fitting basis, which must be specified with the name "rimp2 basis" (see [[Basis]]). For instance,
 +  
   +  
  basis "rimp2 basis"
 +  
  O s; 10000.0 1
 +  
  O s; 1000.0 1
 +  
  O s; 100.0 1
 +  
  ...
 +  
  end
 +  
   +  
  Alternatively, using a standard capability of basis sets ([[Basis]]) another named basis may be associated with the fitting basis. For instance, the following input specifies a basis with the name "small fitting basis" and then defines this to be the "rimp2 basis".
 +  
   +  
  basis "small fitting basis"
 +  
  H s; 10 1
 +  
  H s; 3 1
 +  
  H s; 1 1
 +  
  H s; 0.1 1
 +  
  H s; 0.01 1
 +  
  end
 +  
   +  
  set "rimp2 basis" "small fitting basis"
 +  
   +  
  ==FILE3C  RIMP2 3center integral filename==
 +  
   +  
  The default name for the file used to store the transformed 3center integrals is "$file_prefix$.mo3cint" in the scratch directory. This may be overridden using the FILE3C directive. For instance, to specify the file /scratch/h2o.3c, use this directive
 +  
   +  
  file3c /scratch/h2o.3c
 +  
   +  
  ==RIAPPROX  RIMP2 Approximation==
 +  
   +  
  The type of RI approximation used in the RIMP2 calculation is controlled by means of the RIAPPROX directive. The two possible values are V and SVS (case sensitive), which correspond to the approximations with the same names described in O. Vahtras, J Almlöf, and M. W. Feyereisen, Chem. Phys. Lett. 213, 514518 (1993). The default is V.
 +  
   +  
  ==Advanced options for RIMP2==
 +  
   +  
  These options, which functioned at the time of writing, are not currently supported.
 +  
   +  
  ===Control of linear dependence===
 +  
   +  
  Construction of the RI fit requires the inversion of a matrix of fitting basis integrals which is carried out via diagonalization. If the fitting basis includes near linear dependencies, there will be small eigenvalues which can ultimately lead to nonphysical RIMP2 correlation energies. Eigenvectors of the fitting matrix are discarded if the corresponding eigenvalue is less than <math>mineval</math> which defaults to <math>10^{8}</math>. This parameter may be changed by setting the a parameter in the database. For instance, to set it to <math>10^{10}</math>
 +  
   +  
  set "mp2:fit min eval" 1e10
 +  
   +  
  ===Reference Spin Mapping for RIMP2 Calculations===
 +  
   +  
  The user has the option of specifying that the RIMP2 calculations are to be done with variations of the SCF reference wavefunction. This is accomplished with a SET directive of the form,
 +  
   +  
  set "mp2:reference spin mapping" <integer array default 0>
 +  
   +  
  Each element specified for array is the SCF spin case to be used for the corresponding spin case of the correlated calculation. The number of elements set determines the overall type of correlated calculation to be performed. The default is to use the unadulterated SCF reference wavefunction.
 +  
   +  
  For example, to perform a spinunrestricted calculation (two elements) using the alpha spin orbitals (spin case 1) from the reference for both of the correlated reference spin cases, the SET directive would be as follows,
 +  
   +  
  set "mp2:reference spin mapping" 1 1
 +  
   +  
  The SCF calculation to produce the reference wavefunction could be either RHF or UHF in this case.
 +  
   +  
  The SET directive for a similar case, but this time using the betaspin SCF orbitals for both correlated spin cases, is as follows,
 +  
   +  
  set "mp2:reference spin mapping" 2 2
 +  
   +  
  The SCF reference calculation must be UHF in this case.
 +  
   +  
  The SET directive for a spinrestricted calculation (one element) from the betaspin SCF orbitals using this option is as follows,
 +  
   +  
  set "mp2:reference spin mapping" 2
 +  
   +  
  The SET directive for a spinunrestricted calculation with the spins flipped from the original SCF reference wavefunction is as follows,
 +  
   +  
  set "mp2:reference spin mapping" 2 1
 +  
   +  
  ===Batch Sizes for the RIMP2 Calculation===
 +  
   +  
  The user can control the size of each batch in the transformation and energy evaluation in the MP2 calculation, and consequently the memory requirements and number of passes required. This is done using two SET directives of the following form,
 +  
   +  
  set "mp2:transformation batch size" <integer size default 1>
 +  
  set "mp2:energy batch size" <integer isize jsize default 1 1>
 +  
   +  
  The default is for the code to determine the batch size based on the available memory. Should there be problems with the programdetermined batch sizes, these variables allow the user to override them. The program will always use the smaller of the user's value of these entries and the internally computed batch size.
 +  
   +  
  The transformation batch size computed in the code is the number of occupied orbitals in the ''(occ virfit)'' threecenter integrals to be produced at a time. If this entry is less than the number of occupied orbitals in the system, the transformation will require multiple passes through the twoelectron integrals. The memory requirements of this stage are two global arrays of dimension ''<batch size> x vir x fit'' with the "fit" dimension distributed across all processors (on shellblock boundaries). The compromise here is memory space versus multiple integral evaluations.
 +  
   +  
  The energy evaluation batch sizes are computed in the code from the number of occupied orbitals in the two sets of threecenter integrals to be multiplied together to produce a matrix of approximate fourcenter integrals. Two blocks of integrals of dimension ''(<batch isize> x vir)'' and ''(<batch jsize> x vir)'' by fit are read in from disk and multiplied together to produce ''<batch isize> <batch jsize> vir^2'' approximate integrals. The compromise here is performance of the distributed matrix multiplication (which requires large matrices) versus memory space.
 +  
   +  
  ===Energy Memory Allocation Mode: RIMP2 Calculation===
 +  
   +  
  The user must choose a strategy for the memory allocation in the energy evaluation phase of the RIMP2 calculation, either by minimizing the amount of I/O, or minimizing the amount of computation. This can be accomplished using a SET directive of the form,
 +  
   +  
  set "mp2:energy mem minimize" <string mem_opt default I>
 +  
   +  
  A value of I entered for the string mem_opt means that a strategy to minimize I/O will be employed. A value of C tells the code to use a strategy that minimizes computation.
 +  
   +  
  When the option to minimize I/O is selected, the block sizes are made as large as possible so that the total number of passes through the integral files is as small as possible. When the option to minimize computation is selected, the blocks are chosen as close to square as possible so that permutational symmetry in the energy evaluation can be used most effectively.
 +  
   +  
  ===Local Memory Usage in ThreeCenter Transformation===
 +  
   +  
  For most applications, the code will be able to size the blocks without help from the user. Therefore, it is unlikely that users will have any reason to specify values for these entries except when doing very particular performance measurements.
 +  
   +  
  The size of xf3ci:AO 1 batch size is the most important of the three, in terms of the effect on performance.
 +  
   +  
  Local memory usage in the first two steps of the transformation is controlled in the RIMP2 calculation using the following SET directives,
 +  
   +  
  set "xf3ci:AO 1 batch size" <integer max>
 +  
  set "xf3ci:AO 2 batch size" <integer max>
 +  
  set "xf3ci:fit batch size" <integer max>
 +  
   +  
  The size of the local arrays determines the sizes of the two matrix multiplications. These entries set limits on the size of blocks to be used in each index. The listing above is in order of importance of the parameters to performance, with xf3ci:AO 1 batch size being most important.
 +  
   +  
  Note that these entries are only upper bounds and that the program will size the blocks according to what it determines as the best usage of the available local memory. The absolute maximum for a block size is the number of functions in the AO basis, or the number of fitting basis functions on a node. The absolute minimum value for block size is the size of the largest shell in the appropriate basis. Batch size entries specified for max that are larger than these limits are automatically reset to an appropriate value.
 +  
   +  
  ==Oneelectron properties and natural orbitals==
 +  
   +  
  If an MP2 energy gradient is computed, all contributions are available to form the MP2 linearresponse density. This is the density that when contracted with any spinfree, oneelectron operator yields the associated property defined as the derivative of the energy. Thus, the reported MP2 dipole moment is the derivative of the energy w.r.t. an external magnetic field and is not the expectation value of the operator over the wavefunction. Only dipole moments are printed by the MP2 gradient code, but natural orbitals are produced and stored in the permanent directory with a file extension of ".mp2nos". These may be fed into the [[Propertiesproperty package]] to compute more general properties. Note that the MP2 linear response density matrix is not necessarily positive definite so it is not unusual to see a few small negative natural orbital occupation numbers.
 +  
   +  
  ==References==
 +  
  <references/>
 +  