653. Energy decomposition analysis the manual/multiwfn way -- nwchem

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 

649. N/EDA in GAMESS. 1. Recompiling GAMESS US with NBO6

I need to do energy decomposition analysis (EDA), but only have licenses for Gaussian and NBO6 (i.e. not ADF, turbomole, QChem etc.). NEDA isn't supported by NBO6 with gaussian (afaik).

NWChem, my usual gaussian alternative, doesn't support NBO6 beyond writing a .47 file.

Enter GAMESS US. I've been trying out gamess every few years, but I've found that it's slow and unreliable (very difficult to get SCF convergence) for the systems I work with (polyanions). Some of this may obviously be down to my lack of familiarity with the code -- there are probably plenty of satisified users of GAMESS US.

Either way, NBO6 suppports NEDA with GAMESS US. Also, GAMESS US does two different types  of   EDA: the common Morokuma-Kitaura (MOROKM) one (although only with HF) and an alternative approach by Su and Li that's referred to by GAMESS as LMOEDA (or CMOEDA).

MOROKM and LMOEDA as supported out of the box by GAMESS, but to get it to do NEDA you need to re-link it against NBO. Luckily it's even easier than the instructions in the NBO gamess file (i.e. no need to edit code).

NOTE: I could only link with gfortran 4.9 (jessie). gfortran 6.3 (stretch) failed to link (messages re -fPIC; recompiling gamess with -fPIC didn't solve it).

To compile gamess, see e.g. http://verahill.blogspot.com/2013/06/4xx-gamess-us-2013-r1-on-debian-wheezy.html

Once you've done the ddi/comp and compall steps, edit lked and search for NBO. Change to

set NBO=true
set NBOLIB="/opt/nbo6/bin/gmsnbo.i8.a"

assuming that this location is correct.

Then do lked as in the post above.