Many-body expansion
Compilation
To run MBE
for large systems, the default maximum number of atoms must be increased. This can be done via configure. The maximum number that can be entered here is 1000. To increase this further the configure script must be modified. The default number of records should also be increased. This should typically be about 4 times the number of monomers in the system, as multipoles, polarizabilities and density-fitting information are held separate records for each monomer.
Units
All lengths in the input file should be given in Angstroms.
Incremental Monte-Carlo
LMAX_M
Angular momentum of multipoles to be usedLMAX_P
Angular momentum of polarizabilities to be usedPRINT
Level of information to outputIMCDUMP
At end of a calculation dump various levels of information: $>=1$ dumps properties and one-body energies, $>=2$ dumps all two-body energies, $>=3$ all three-body energiesMAXLEV
Maximum level of MBE. Default is to include dimers (=2). For IMC maximum is trimers. For a manybody analysis can go up to 5.CUTOFF
Distance within which to do QM calculations. For 3 body calculations, all dimer separations must be within this distance.THRDENOV
Distance within which to calculate the overlap for damping. By default the same asCUTOFF
.MC
Do a Monte-Carlo simulation if 1 ( Default 0)MCMAX
Maximum number of MC steps to performMCTEMP
Temperature of a MC simulationNPT
Do a NPT simulationPRESSURE
Pressure for a NPT MC simulationPRFREQ
How often to print output during a MC simulationRESFREQ
How often to save restart information (may not work!)ENG_PROC
Name of energy procedure for all levels (if the same is to be used for all)ENG_PROC1
Name of energy procedure for monomersENG_PROC2
Name of energy procedure for dimersENG_PROC3
Name of energy procedure for trimersPROP_PROC
Name of property procedureRESTART
If we are doing a restart (may not work)DAMPING
Use damping (Tang-Toennies with empirical factor of 1.94 by default)DAMPFAC
Empirical factor to use in damping functionDF_SET
Density-fitting set to use in damping. Density damping is only switched on if this is specifiedDF_CON
Optional specification ofCONTEXT
forDF_SET
, e.g.JKFIT
orMP2FIT
.JKFIT
by defaultDFMAXL
Optional specification of maximum angular momentum functions to use inDF_SET
DENREC
Location where density has been stored for calculating distributed multipolesGEOM_LOC
Location where geometry output from MC should be storedANALYSE
Only do analysis and no calculationRDF_TYPES
Type of RDFs to calculateRDF_VOL
Volume for non-periodic RDFRDF_NBIN
Number of bins for RDFRDF_DBIN
Width of binRDF_EQ
Number of equilibration steps to ignore when calculating RDFMANYBODY
Just do a many-body analysis. =1 does normal many-body analysis and =2 does full counterpoise corrected.SCALTYPE
Scale box before starting calculation. =A for side of box in Angstroms, =B for side of box in Bohr, =N for number density, =D for density in kg/m3SCALE
Size of scaled box in whatever type has been specified inSCALTYPE
DTRANS
Initial size of translational move in AngstromDROT
Initial size of rotational move in degreesDSTR
Initial size of bond stretch in AngstromDBEND
Initial size of bond angle bend in degreesDBOX
Initial size of box move in Angstroms (for NPT simulations)DINTRA
Bias for intramolecular moves (this plusDRIGID
will be equal to one.DINTRA
takes priority, if both specified)DRIGID
Bias for rigid-body moves (this plusDINTRA
will be equal to one.DINTRA
takes priority, if both specified)DVOL
Bias for volume moves versus molecular moves in NPT simulationsTARGINT
Acceptance/rejection ratio target for intramolecular movesTARGRIG
Acceptance/rejection ratio target for rigid-body movesTARGVOL
Acceptance/rejection ratio target for volume moves in NPT simulationsMONFILE
File to specify connectivity of system which overrides automatic connectivity subroutine. Should not be used with MC simulations, but OK for energy.CHGFILE
File to specify if any of the monomers are charged.ANISODISP
Are we using anisotropic dispersion integralsDISPFILE
File to specify dispersion coefficients or integrals, bohr for isotropic and anisotropicINDMETH
Method to use for self-consistent induction calculation. =0 don’t iterate, =1 iterate and use Ewald only on first iteration for PBC, =2 iterate and use Ewald at every step for PBC, =3 use Lanczos algorithm which uses Ewald on first iteration only.ITNUM
Starting value for the iteration number.NBODY
Include n-body energy (default =1)CLASSICAL
Include classical energies (default =1)SWITCH
Use a switching function. =0 no switching, =1 quintic splines switching, (=2 GAP switching)PROPFILE
File which contains multipoles and polarizabilities when doing a MaxLev=0 calculation, i.e. no QM calculations being doneDFPROC
Name of density-fitting procedure when only model being used and properties being input. Done only once for each type of monomerEXCH_K
Factor by which to multiply the overlap to give the exchange energyMANYBODY
Do a many-body analysis on the system. =1 do normal MBE, =2 do MBE in basis of cluster, i.e. properly counterpoise-correctedDMAMON
Do a distributed multipole analysis on the whole system and also output the multipoles at the centres-of-mass of each monomer, due to contributions from that monomer only. Should be comparable to induced multipoles.CLOSEST
=$N$ Output the closest $N$ molecules (possibly within minimum image convention if PBC being used) to a molecule chosen at random. Useful for creating clustersSUPERCELL
Output a supercell specified as a string ’na:nb:nc’ using the PBC specifications, e.g.SUPERCELL=’1:1:1’
creates an array of 8 of the original cells in a ’cubic’ arrangementSUPERFILE
File to which the supercell geometry is output
EWALD directive
GAMMA
Exponent of screening gaussian in Ewald summationKCUT
Cutoff radius for radius of k-vectorsRCUT
Cutoff radius for real-space ewaldLCUT
Cutoff for angular momentum in multipolar Ewald. Higher-order terms probably convergent in real-space radiusEPS
Permittivity to use for surface term in Ewald
Examples
Obtaining TDHF and UCHF dispersion coefficients
The example below demonstrates how to get T. Korona’s CC-SAPT program to output TDHF and UCHF dispersion coefficients. The UCHF coefficients are the relevant ones for an MP2 calculation. This example gives dispersion coefficients for water-water, ammonia-ammonia and water-ammonia. To use them in a calculation all of the integrals listed in the output for a given interaction should be put into a file in exactly the format they are output, i.e. four indices and an integral. This can then be read by the MBE
program.
- examples/h2o_nh3_disp.inp
memory,100,m if(NPROC_MPP.gt.1) then skipped end if gthresh,energy=1.d-10 ! don't use symmetry - can mess things up symmetry,nosym geometry={ H O 1 R H 2 R 1 A } R=0.966443838 Angstrom A=103.84335748 Degree gexpec,nspmlt3 ileden=0 iledenp=0 ! first call just sets up CCSAPT program dispgg=15 set,CC_NORM_MAX=50 basis=sto-3g hf;maxit,0 {ccsd,check=0 polari,nspmlt3 orbital,ignore_error=1;maxit,3 cprop,ccsapt=3} ! now calculate the dispersion integrals for water basis=avdz hf;save,scfa tdhf=tdhfa tdhfdisp=tdhfdispa {ccsd,check=0 polari,nspmlt3 orbital,ignore_error=1;maxit,3 cprop,ccsapt=5,dispcoef=dispgg,dispomega=0.3} ! calculate dispersion integrals for ammonia and mixed coefficient for NH3-H2O symmetry,nosym geometry={ N X1 N RX H1 N RNH X1 DUMB H2 N RNH X1 DUMB H1 DIHE 0 H3 N RNH X1 DUMB H1 -DIHE 0 } RX= 1.02 Angstrom RNH= 1.02071411 Angstrom DIHE= 120. Degree DUMB= 112.48619220 Degree ileden=0 iledenp=0 gexpec,nspmlt3 dispgg=15 set,CC_NORM_MAX=50 basis=sto-3g hf;maxit,0 {ccsd,check=0 polari,nspmlt3 orbital,ignore_error=1;maxit,3 cprop,ccsapt=3} basis=avdz hf;save,2050.2 tdhf=tdhfa tdhfdisp=tdhfdispa {ccsd,check=0 polari,nspmlt3 orbital,ignore_error=1;maxit,3 cprop,ccsapt=5,dispcoef=dispgg,dispomega=0.3}
Many-body analysis of water hexamer
This example performs a many-body analysis of the water hexamer up to 3-body contributions. The routines go up to a maximum of 5-body contributions. By specifying MANYBODY=2
a full counterpoise corrected many-body analysis would be done.
- examples/water6_mbe.inp
skipped ! bug5219 PROC mbeeng {df-hf;start,atden} {df-lmp2,interact=1 enepart} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geometry={ 18 DF-LMP2/AVDZ ENERGY=-457.63602408 O -2.1233226309 -1.7688858791 0.1524109399 H -2.4367986815 -2.2103390376 0.9513842141 H -1.2045230633 -2.0982518768 0.0330050109 O 2.5947061994 -0.9537736740 0.1422730545 H 2.4193021370 0.0068583194 0.0266664332 H 3.1418238405 -1.0061708640 0.9356521604 O -2.5946186870 0.9534257806 -0.1546242536 H -3.1315700574 1.0016502753 -0.9551345770 H -2.4201260402 -0.0065436988 -0.0322861934 O 2.1219816165 1.7703156035 -0.1450981681 H 1.2027719526 2.0987725437 -0.0263944582 H 2.4388614383 2.2190153349 -0.9386745371 O -0.4725613825 2.7226532886 0.1512286611 H -0.7014891628 3.2160843341 0.9484803611 H -1.2168653397 2.0916787132 0.0293211623 O 0.4703934181 -2.7233099281 -0.1463373783 H 0.6968166194 -3.2213522313 -0.9414479081 H 1.2152177361 -2.0918266148 -0.0304245102 } basis,avdz {mbe,manybody=1,maxlev=3}
Two-body-plus-polarization treatment of water hexamer
This demonstrates a many-body evaluation of the energy of the cluster using the long-range model for well-separated dimers and the polarization model for many-body (beyond two-body) effects. Two procedures are defined, one for the energy and one for the properties. Properties up to quadrupoles (LMAX_M=2, LMAX_P=2
) are used. MBEENG
and MBEPROP
are the default names for these procedures and hence do not need to be specified on the input line. MAXLEV=2
indicates that up to dimers should be calculated using MBEENG
and within CUTOFF=4.5
(in Angstroms). Isotropic dispersion coefficients of Wormer et al. are used by default for water.
- examples/water6_mbe2.inp
skipped ! bug5219 PROC mbeeng !default procedure name {hf} ENDPROC PROC mbeprop !default procedure name {hf polarizability,nspmlt2} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geometry={ 18 Water hexamer O -2.1233226309 -1.7688858791 0.1524109399 H -2.4367986815 -2.2103390376 0.9513842141 H -1.2045230633 -2.0982518768 0.0330050109 O 2.5947061994 -0.9537736740 0.1422730545 H 2.4193021370 0.0068583194 0.0266664332 H 3.1418238405 -1.0061708640 0.9356521604 O -2.5946186870 0.9534257806 -0.1546242536 H -3.1315700574 1.0016502753 -0.9551345770 H -2.4201260402 -0.0065436988 -0.0322861934 O 2.1219816165 1.7703156035 -0.1450981681 H 1.2027719526 2.0987725437 -0.0263944582 H 2.4388614383 2.2190153349 -0.9386745371 O -0.4725613825 2.7226532886 0.1512286611 H -0.7014891628 3.2160843341 0.9484803611 H -1.2168653397 2.0916787132 0.0293211623 O 0.4703934181 -2.7233099281 -0.1463373783 H 0.6968166194 -3.2213522313 -0.9414479081 H 1.2152177361 -2.0918266148 -0.0304245102 } basis,avdz {mbe,maxlev=2,cutoff=4.5,lmax_m=2,lmax_p=2}
Three-body-plus-polarization treatment of water hexamer
Same as previous but with MAXLEV=3
so that trimers are also calculated using MBEENG
. Note that SWITCH=0
has been used as no 3-body switching function has been implemented and ‘unswitched’ 3-body energies should not be mixed with ‘switched’ 2-body ones.
- examples/water6_mbe3.inp
skipped ! bug5219 PROC mbeeng !default procedure name {hf} ENDPROC PROC mbeprop !default procedure name {hf polarizability,nspmlt2} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geometry={ 18 Water hexamer O -2.1233226309 -1.7688858791 0.1524109399 H -2.4367986815 -2.2103390376 0.9513842141 H -1.2045230633 -2.0982518768 0.0330050109 O 2.5947061994 -0.9537736740 0.1422730545 H 2.4193021370 0.0068583194 0.0266664332 H 3.1418238405 -1.0061708640 0.9356521604 O -2.5946186870 0.9534257806 -0.1546242536 H -3.1315700574 1.0016502753 -0.9551345770 H -2.4201260402 -0.0065436988 -0.0322861934 O 2.1219816165 1.7703156035 -0.1450981681 H 1.2027719526 2.0987725437 -0.0263944582 H 2.4388614383 2.2190153349 -0.9386745371 O -0.4725613825 2.7226532886 0.1512286611 H -0.7014891628 3.2160843341 0.9484803611 H -1.2168653397 2.0916787132 0.0293211623 O 0.4703934181 -2.7233099281 -0.1463373783 H 0.6968166194 -3.2213522313 -0.9414479081 H 1.2152177361 -2.0918266148 -0.0304245102 } basis,avdz {mbe,maxlev=3,cutoff=4.5,lmax_m=2,lmax_p=2,switch=0}
Specifying isotropic dispersion coefficients
Similar to the example in Subsection Two-body-plus-polarization treatment of water hexamer, but here we are specifying the isotropic dispersion coefficients which should be used. The DISPFILE
variable specifies the file where these should be found. It should be a list of the relevant interactions and the $C_6$, $C_8$ and $C_{10}$ coefficients. For example
''%%OH2:OH2 46.0 800.0 10000.0 %%''\\
where molecules should be listed by elements of decreasing atomic mass. If a charged species is being used, the charge should be included in parentheses, e.g. F(-1)
. The file may contain interactions not relevant to the system under consideration, so a library of dispersion coefficients may be built up. The input file is
- examples/water6_iso.inp
skipped ! bug5219 PROC mbeeng !default procedure name {hf} ENDPROC PROC mbeprop !default procedure name {hf polarizability,nspmlt2} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geometry={ 18 Water hexamer O -2.1233226309 -1.7688858791 0.1524109399 H -2.4367986815 -2.2103390376 0.9513842141 H -1.2045230633 -2.0982518768 0.0330050109 O 2.5947061994 -0.9537736740 0.1422730545 H 2.4193021370 0.0068583194 0.0266664332 H 3.1418238405 -1.0061708640 0.9356521604 O -2.5946186870 0.9534257806 -0.1546242536 H -3.1315700574 1.0016502753 -0.9551345770 H -2.4201260402 -0.0065436988 -0.0322861934 O 2.1219816165 1.7703156035 -0.1450981681 H 1.2027719526 2.0987725437 -0.0263944582 H 2.4388614383 2.2190153349 -0.9386745371 O -0.4725613825 2.7226532886 0.1512286611 H -0.7014891628 3.2160843341 0.9484803611 H -1.2168653397 2.0916787132 0.0293211623 O 0.4703934181 -2.7233099281 -0.1463373783 H 0.6968166194 -3.2213522313 -0.9414479081 H 1.2152177361 -2.0918266148 -0.0304245102 } basis,avdz {mbe,maxlev=2,cutoff=4.5,mc=0,lmax_m=2,lmax_p=2,anisodisp=0,dispfile='disp_iso.dat'}
and the required dispersion file is:
- examples/disp_iso.dat
OH2:OH2 40.196 496.889 8263.5
Specifying anisotropic dispersion coefficients
Again this is controlled by the DISPFILE
variable, but the file now has a different format. For each type of interaction three filenames should be specified. The first should contain the dispersion integrals and the second two the reference geometries at which they were calculated, e.g.
OH2:OH2 h2oh2odisp.dat h2o.xyz h2o.xyz
The geometries should be given in standard XYZ format. The input file is:
- examples/water6_aniso.inp
skipped ! bug5219 PROC mbeeng !default procedure name {hf} ENDPROC PROC mbeprop !default procedure name {hf polarizability,nspmlt2} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geometry={ 18 Water hexamer O -2.1233226309 -1.7688858791 0.1524109399 H -2.4367986815 -2.2103390376 0.9513842141 H -1.2045230633 -2.0982518768 0.0330050109 O 2.5947061994 -0.9537736740 0.1422730545 H 2.4193021370 0.0068583194 0.0266664332 H 3.1418238405 -1.0061708640 0.9356521604 O -2.5946186870 0.9534257806 -0.1546242536 H -3.1315700574 1.0016502753 -0.9551345770 H -2.4201260402 -0.0065436988 -0.0322861934 O 2.1219816165 1.7703156035 -0.1450981681 H 1.2027719526 2.0987725437 -0.0263944582 H 2.4388614383 2.2190153349 -0.9386745371 O -0.4725613825 2.7226532886 0.1512286611 H -0.7014891628 3.2160843341 0.9484803611 H -1.2168653397 2.0916787132 0.0293211623 O 0.4703934181 -2.7233099281 -0.1463373783 H 0.6968166194 -3.2213522313 -0.9414479081 H 1.2152177361 -2.0918266148 -0.0304245102 } basis,avdz {mbe,maxlev=2,cutoff=4.5,mc=0,lmax_m=2,lmax_p=2,anisodisp=1,dispfile='disp.dat'}
The additional files that are required are: disp.dat, h2o_tip.xyz and h2o_tip_uchf.disp.
Using MP2 properties
As the example in Subsection Two-body-plus-polarization treatment of water hexamer but using MP2 properties. We have now also specified different procedures for monomer and dimer energies. This is important for correlated energies as the specifying the correct number of core orbitals is important. The examples require auxiliary files disp.dat, h2o_tip.xyz and h2o_tip_uchf.disp.
- examples/water6_mp2.inp
skipped ! bug5219 PROC mbeeng1 core,1 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC mbeeng2 core,2 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC mbeprop core,0 {hf;start,atden polarizability,NSPMLT2} {mp2;dm,11000.2} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geometry={ 18 Water hexamer O -2.1233226309 -1.7688858791 0.1524109399 H -2.4367986815 -2.2103390376 0.9513842141 H -1.2045230633 -2.0982518768 0.0330050109 O 2.5947061994 -0.9537736740 0.1422730545 H 2.4193021370 0.0068583194 0.0266664332 H 3.1418238405 -1.0061708640 0.9356521604 O -2.5946186870 0.9534257806 -0.1546242536 H -3.1315700574 1.0016502753 -0.9551345770 H -2.4201260402 -0.0065436988 -0.0322861934 O 2.1219816165 1.7703156035 -0.1450981681 H 1.2027719526 2.0987725437 -0.0263944582 H 2.4388614383 2.2190153349 -0.9386745371 O -0.4725613825 2.7226532886 0.1512286611 H -0.7014891628 3.2160843341 0.9484803611 H -1.2168653397 2.0916787132 0.0293211623 O 0.4703934181 -2.7233099281 -0.1463373783 H 0.6968166194 -3.2213522313 -0.9414479081 H 1.2152177361 -2.0918266148 -0.0304245102 } basis,avdz {mbe,maxlev=2,cutoff=4.5,mc=0,lmax_m=2,lmax_p=2,eng_proc1='MBEENG1',eng_proc2='MBEENG2',anisodisp=1,dispfile='disp.dat'}
A second example is provided which produces more verbose output.
- examples/water6_mp2_moreinfo.inp
skipped ! bug5219 PROC MBEENG1 core,1 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC MBEENG2 core,2 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC MBEPROP core,0 {hf;start,atden polarizability,NSPMLT2} {mp2;dm,11000.2} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geometry={ 18 Water hexamer O -2.1233226309 -1.7688858791 0.1524109399 H -2.4367986815 -2.2103390376 0.9513842141 H -1.2045230633 -2.0982518768 0.0330050109 O 2.5947061994 -0.9537736740 0.1422730545 H 2.4193021370 0.0068583194 0.0266664332 H 3.1418238405 -1.0061708640 0.9356521604 O -2.5946186870 0.9534257806 -0.1546242536 H -3.1315700574 1.0016502753 -0.9551345770 H -2.4201260402 -0.0065436988 -0.0322861934 O 2.1219816165 1.7703156035 -0.1450981681 H 1.2027719526 2.0987725437 -0.0263944582 H 2.4388614383 2.2190153349 -0.9386745371 O -0.4725613825 2.7226532886 0.1512286611 H -0.7014891628 3.2160843341 0.9484803611 H -1.2168653397 2.0916787132 0.0293211623 O 0.4703934181 -2.7233099281 -0.1463373783 H 0.6968166194 -3.2213522313 -0.9414479081 H 1.2152177361 -2.0918266148 -0.0304245102 } basis,avdz {mbe,maxlev=2,cutoff=4.5,mc=0,lmax_m=2,lmax_p=2,eng_proc1='MBEENG1',eng_proc2='MBEENG2',anisodisp=1,dispfile='disp.dat',print=3,imcdump=2}
Mixed water–ammonia cluster
A mixed cluster of 10 water molecules and two ammonia molecules. There is now anisotropic dispersion information listed for multiple interactions. The dispersion file disp_h2o_nh3.dat specifies monomer reference geometries (h2o.xyz and nh3.xyz) and UCHF dispersion coefficients in h2o_nh3_uchf.dat, h2o_uchf.dat and nh3_uchf.dat. The geometry is specified in water10_ammonia2.xyz and the input file is below.
- examples/water10_ammonia2.inp
skipped ! bug5219 memory,40,m PROC MBEENG1 core,1 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC MBEENG2 core,2 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC MBEPROP core,0 {hf;start,atden polarizability,nspmlt3} {mp2;dm,11000.2} ENDPROC symmetry,nosym geomtyp=xyz geom=water10_ammonia2.xyz basis,avdz {mbe,maxlev=2,cutoff=4.5,lmax_m=3,lmax_p=3,damping=1,df_set='CC-PVTZ(D/P)',dispfile='disp.dat',eng_proc1='MBEENG1',eng_proc2='MBEENG2',switch=0}
Single calculation on 64 water molecules in periodic boundary conditions
The PBC
command specifies that periodic boundary conditions should be invoked. The ewald
directive is also present so that periodic electrostatic and induction energies are included. Without this directive everything would simply be done in the minimum image convention. Damping is used and since a fitting set has been specified using DF_SET
, damping will be done using density overlap of monomers. The required auxiliary files for this example are:
- geometry: water64.xyz
The input file is below.
- examples/water64_pbc.inp
skipped ! bug5219 if(NPROC_MPP.gt.1) then skipped end if ! define the procedure to be run on the monomers PROC MBEENG1 core,1 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC MBEENG2 core,2 {df-hf;start,atden} {df-lmp2,interact=1 pipek,delete=1 enepart} ENDPROC PROC MBEPROP core,1 {hf;start,atden polarizability,nspmlt2} {mp2;dm,11000.2} ENDPROC ! input ensemble geometry symmetry,nosym orient,noorient geomtyp=xyz geom=water64.xyz basis=avdz {pbc,latt_type='CUBIC',boxsize=12.4297299} {mbe,maxlev=2,print=1,lmax_m=2,lmax_p=2,cutoff=4.5,damping=1,DF_SET='CC-PVTZ(D/P)',eng_proc1='MBEENG1',eng_proc2='MBEENG2',dispfile='disp.dat',anisodisp=1 ewald}