Multiwfn forum

b322qr

Member

Registered: 2024-04-30

Posts: 7

sobEDA Unstable Wavefunction

I am using sobEDA on an actinide system and after it has completed, I run a stability check in Gaussian using the "final.chk" checkpoint file, which then shows that the calculation has an unstable wavefunction. What should I do in order to get the correct energy decomposition analysis with the stable wavefunction? Would I be able to run another second single point calculation using the stable wavefunction as a guess wavefunction and extract the EDA parameters from there?

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 1,830

Website

Re: sobEDA Unstable Wavefunction

This situation is very special. I can only recommend to hack the sobEDA.sh. Please read through the script to understand its logic, and then properly modify it. The script is not long and not complex. If you cannot understand some codes in the script, you are welcomed to post here. In fact, you can fully manually perform each step of the script.

b322qr

Member

Registered: 2024-04-30

Posts: 7

Re: sobEDA Unstable Wavefunction

Thank you, professor Lu. I appreciate your offer for help. I'm grateful that your code already creates the checkpoint file for each step of the calculation, making this process much easier. I would like your insight on my logic as I try to understand how to use the energy from the stable wavefunction to extract the information necessary for EDA. Most importantly, the "final" file.

I see at the bottom of the sobEDA.sh script that the only energy values that depend on the final energy are the total (dE_tot), dispersion (dE_disp), Pauli repulsion energy (dE_rep), and orbital energy (dE_orb). For reference, these are the lines of code I am looking at in the sobEDA.sh script. 

##### Show summary of interaction #####
dE_tot=$E_tot_final
dE_els=$E_els_promol
dE_x=$E_x_promol
dE_c=$E_c_promol
dE_disp=$E_disp_final
for ((i=1;i<=$nfrag;i=i+1))
do
  dE_tot=`echo "$dE_tot-(${E_tot[$i]})" | bc -l`
  dE_els=`echo "$dE_els-(${E_els[$i]})" | bc -l`
  dE_x=`echo "$dE_x-(${E_x[$i]})" | bc -l`
  dE_c=`echo "$dE_c-(${E_c[$i]})" | bc -l`
  dE_disp=`echo "$dE_disp-(${E_disp[$i]})" | bc -l`
done
dE_orb=`echo "$E_tot_final-($E_tot_frz)" | bc -l`
dE_rep=`echo "$E_tot_frz-($E_tot_promol)" | bc -l`

The final total energy (E_tot_final) looks for the energy in the line that says "SCF Done" in final.out, dispersion energy (E_disp_final) looks for the energy in the line that says "Dispersion energy" in final.out, frozen energy (E_tot_frz) looks for the energy shown in the first cycle in final.out (which final.gjf reads from the promol.chk file), and E_tot_promol looks for the energy in the line that says "SCF Done" in promol.out.

After running a single point calculation (let's call this new single point calculation "final2") using the same route section as final.gjf, but reading the guess from the stable wavefunction, I receive an output that finishes after only 1 cycle (I guess because it is reading from an already optimized, stable wavefunction). With the knowledge of the way the parameters of EDA are calculated by the program, I compare between final and final2 the dispersion energy (it does not change, and therefore dE_disp would not change) and the final energy printed after "SCF Done" (it becomes more negative because of the stable wavefunction).

If I try to perform the EDA on final2.out, obviously I will get erroneous results because the frozen energy is the same as the final energy (it finishes in one cycle). So, my next thought is to modify the line in the original final.out that says "SCF Done" and replace it with the energy from that same line in final2.out. (i.e. change  "SCF Done:  E(UPBE1PBE) =  XXX.XXX" to " SCF Done:  E(UPBE1PBE) =  YYY.YYY") That way the final energy represents the energy from the stable wavefunction. Although, this is assuming the frozen energy would remain the same.

Please let me know if my logic is correct and if it is reasonable to replace the energy value shown after "SCF Done" in final.out with the energy from final2.out.

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 1,830

Website

Re: sobEDA Unstable Wavefunction

You can see the following energy terms are extracted from final.out

E_tot_final=`grep "SCF Done" final.out | awk -F = '{print $2}'| awk '{print $1}'`
E_T_final=`grep "ET=" final.out | awk -F = '{print $2}' | awk '{print $1}'`
E_x_final=`grep "ENTVJ=" final.out | awk -F = '{print $3}' | awk '{print $1}'`
E_c_final=`grep "ENTVJ=" final.out | awk -F = '{print $4}' | awk '{print $1}'`
...
E_disp_final=`grep "Dispersion energy" final.out | awk -F = '{print $2}'| awk '{print $1}'`
...
E_tot_frz=`grep " E=" final.out | head -1 | awk '{print $2}'`

Before running these and subsequent lines in the sobEDA.sh, you only need to replace the values after "SCF Done", "ET=", "ENTVJ=" in final.out with those in final2.out.