I have a very large system (390 atoms, 3918 functions, 6474 primitives) where I want to analyse the bonding. Whereas I can reduce the size of the system a little bit, there's a large conjugated ad charged system in the middle which I can't really reduce. Either way, when I use GAMESS US to do NEDA, the calc seems to hang for days without ever progressing, and LMOEDA/CMOEDA keep running out of memory.
I recently had a look at Multiwfn, and section 4.100.8 in the
manual shows how to do simple EDA as a multistep computation. The example uses multiwfn to input initial fragment wavefunctions to compute the DE_orb with Gaussian. Incidentally, this is something which is very easy to do with nwchem without using multiwfn.
I'll use NH3..BH3 as the example at RHF/6-31G*.
Nwchem:
1. Optimise NH3..BH3
scratch_dir /home/me/scratch
Title "NH3BH3-nw"
Start NH3BH3-nw
charge 0
geometry noautosym noautoz units angstrom
N 0.0720500 -0.00961700 -0.336156
H 0.871540 0.292859 -0.862886
H -0.685935 0.618297 -0.534511
H -0.187686 -0.922713 -0.663436
B 0.415920 -0.0366410 1.31774
H 0.709693 1.10584 1.58004
H -0.612009 -0.411733 1.83072
H 1.33214 -0.818018 1.41958
end
basis "ao basis" cartesian print
B library "6-31G*"
H library "6-31G*"
N library "6-31G*"
END
scf
RHF
nopen 0
end
task scf optimize
Energy=-82.61181818
2. Run SE calcs on the BH3 and NH3 fragments:
scratch_dir /home/me/scratch
Title "BH3-nw"
Start BH3-nw
charge 0
geometry noautosym noautoz units angstrom
B 0.192902 -0.0151808 0.928551
H 0.486935 1.12724 1.19093
H -0.834959 -0.390362 1.44154
H 1.10930 -0.796437 1.03042
end
basis "ao basis" cartesian print
B library "6-31G*"
H library "6-31G*"
END
scf
RHF
nopen 0
vectors output bh3.movecs
end
task scf energy
Energy=-26.368337779376
and
scratch_dir /home/me/scratch
Title "NH3-nw"
Start NH3-nw
charge 0
geometry noautosym noautoz units angstrom
N -0.150737 0.0117141 -0.725185
H 0.648571 0.314333 -1.25225
H -0.908696 0.639770 -0.923690
H -0.410582 -0.901249 -1.05260
end
basis "ao basis" cartesian print
N library "6-31G*"
H library "6-31G*"
END
scf
RHF
nopen 0
vectors output nh3.movecs
end
task scf energy
Energy=-56.184296916045
3. Finally, use the two movecs created in step 2:
scratch_dir /home/andy/scratch
Title "NH3BH3-nw"
Start NH3BH3-nw
charge 0
geometry noautosym noautoz units angstrom
N 0.0720500 -0.00961700 -0.336156
H 0.871540 0.292859 -0.862886
H -0.685935 0.618297 -0.534511
H -0.187686 -0.922713 -0.663436
B 0.415920 -0.0366410 1.31774
H 0.709693 1.10584 1.58004
H -0.612009 -0.411733 1.83072
H 1.33214 -0.818018 1.41958
end
basis "ao basis" cartesian print
B library "6-31G*"
H library "6-31G*"
N library "6-31G*"
END
scf
RHF
nopen 0
vectors fragment nh3.movecs bh3.movecs output nh3bh3.movecs
end
task scf
4. Parse the output from step 4:
iter energy gnorm gmax time
----- ------------------- --------- --------- --------
1 -82.5357150919 7.36D-01 2.88D-01 0.1
2 -82.6078664771 2.30D-01 5.23D-02 0.1
3 -82.6117699706 2.03D-02 7.47D-03 0.1
4 -82.6118181287 2.23D-04 5.79D-05 0.1
5 -82.6118181326 2.51D-06 7.28D-07 0.1
Final RHF results
------------------
Total SCF energy = -82.611818132574
One-electron energy = -190.292457149391
Two-electron energy = 67.248334359392
Nuclear repulsion energy = 40.432304657425
Time for solution = 0.1s
So, according to the Multiwfn
Manual at 4.100.8, using the values from above:
DEtot=-82.61181818-(-26.368337779376-56.184296916045)= -37 kcal/mol (-155 kJ/mol)
DEorb=-82.611818132574-(-82.5357150919)= -48 kcal/mol (-200 kJ/mol)
DEsteric=DEtot-DEorb= 11 kcal/mol (45 kJ/mol)
This is essentially the Kitaura-Morokuma method.
See e.g. Frenking et al.in
Energy Decomposition Analysis on page 44. Eq 2 defines Eint in the same way DEtot is defined above, and Eq 7 is the same Eorb as here.
DEstruc here is then DEelstat + DEPauli.
How to resolve these two factors from one another, is a problem for another day.
You can also run the calcs using a single input file:
scratch_dir /home/me/scratch
Title "NH3BH3-nw-eda"
Start NH3BH3-nw-eda
echo
charge 0
geometry molecule noautosym noautoz units angstrom
N -0.150737 0.0117141 -0.725185
H 0.648571 0.314333 -1.25225
H -0.908696 0.639770 -0.923690
H -0.410582 -0.901249 -1.05260
B 0.192902 -0.0151808 0.928551
H 0.486935 1.12724 1.19093
H -0.834959 -0.390362 1.44154
H 1.10930 -0.796437 1.03042
end
geometry ammonia noautosym noautoz units angstrom
N -0.150737 0.0117141 -0.725185
H 0.648571 0.314333 -1.25225
H -0.908696 0.639770 -0.923690
H -0.410582 -0.901249 -1.05260
end
geometry borohydride noautosym noautoz units angstrom
B 0.192902 -0.0151808 0.928551
H 0.486935 1.12724 1.19093
H -0.834959 -0.390362 1.44154
H 1.10930 -0.796437 1.03042
end
basis "ao basis" cartesian print
N library "6-31G*"
B library "6-31G*"
H library "6-31G*"
END
set geometry ammonia
scf
vectors output ammonia.movecs
end
task scf
set geometry borohydride
scf
vectors output borohydride.movecs
end
task scf
set geometry molecule
scf
vectors input fragment ammonia.movecs borohydride.movecs output molecule.movecs
end
task scf