Many-body expansion

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. To compile for parallel execution the MPPX parallel mode should be chosen.

All lengths in the input file should be given in Angstroms.

  • LMAX_M Angular momentum of multipoles to be used
  • LMAX_P Angular momentum of polarizabilities to be used
  • PRINT Level of information to output
  • IMCDUMP 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 energies
  • MAXLEV 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 as CUTOFF.
  • MC Do a Monte-Carlo simulation if 1 ( Default 0)
  • MCMAX Maximum number of MC steps to perform
  • MCTEMP Temperature of a MC simulation
  • NPT Do a NPT simulation
  • PRESSURE Pressure for a NPT MC simulation
  • PRFREQ How often to print output during a MC simulation
  • RESFREQ 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 monomers
  • ENG_PROC2 Name of energy procedure for dimers
  • ENG_PROC3 Name of energy procedure for trimers
  • PROP_PROC Name of property procedure
  • RESTART 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 function
  • DF_SET Density-fitting set to use in damping. Density damping is only switched on if this is specified
  • DF_CON Optional specification of CONTEXT for DF_SET, e.g. JKFIT or MP2FIT. JKFIT by default
  • DFMAXL Optional specification of maximum angular momentum functions to use in DF_SET
  • DENREC Location where density has been stored for calculating distributed multipoles
  • GEOM_LOC Location where geometry output from MC should be stored
  • ANALYSE Only do analysis and no calculation
  • RDF_TYPES Type of RDFs to calculate
  • RDF_VOL Volume for non-periodic RDF
  • RDF_NBIN Number of bins for RDF
  • RDF_DBIN Width of bin
  • RDF_EQ Number of equilibration steps to ignore when calculating RDF
  • MANYBODY 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/m3
  • SCALE Size of scaled box in whatever type has been specified in SCALTYPE
  • DTRANS Initial size of translational move in Angstrom
  • DROT Initial size of rotational move in degrees
  • DSTR Initial size of bond stretch in Angstrom
  • DBEND Initial size of bond angle bend in degrees
  • DBOX Initial size of box move in Angstroms (for NPT simulations)
  • DINTRA Bias for intramolecular moves (this plus DRIGID will be equal to one. DINTRA takes priority, if both specified)
  • DRIGID Bias for rigid-body moves (this plus DINTRA will be equal to one. DINTRA takes priority, if both specified)
  • DVOL Bias for volume moves versus molecular moves in NPT simulations
  • TARGINT Acceptance/rejection ratio target for intramolecular moves
  • TARGRIG Acceptance/rejection ratio target for rigid-body moves
  • TARGVOL Acceptance/rejection ratio target for volume moves in NPT simulations
  • MONFILE 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 integrals
  • DISPFILE File to specify dispersion coefficients or integrals, bohr for isotropic and anisotropic
  • INDMETH 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 done
  • DFPROC Name of density-fitting procedure when only model being used and properties being input. Done only once for each type of monomer
  • EXCH_K Factor by which to multiply the overlap to give the exchange energy
  • MANYBODY Do a many-body analysis on the system. =1 do normal MBE, =2 do MBE in basis of cluster, i.e. properly counterpoise-corrected
  • DMAMON 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 clusters
  • SUPERCELL 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’ arrangement
  • SUPERFILE File to which the supercell geometry is output
  • GAMMA Exponent of screening gaussian in Ewald summation
  • KCUT Cutoff radius for radius of k-vectors
  • RCUT Cutoff radius for real-space ewald
  • LCUT Cutoff for angular momentum in multipolar Ewald. Higher-order terms probably convergent in real-space radius
  • EPS Permittivity to use for surface term in Ewald

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}

There are testjobs already present which show how to calculate dispersion coefficients using CC-SAPT. These should be modified to use the spherical harmonic multipole operators on the EXPEC and POLARI cards to produce the required integrals.

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}

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}

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}

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

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.

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}

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}

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:

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}