Skip to content

Reliable explanation and ways to solve the error messages for the Gaussian quantum chemistry program. This is a per request English translation for a long blog in Chinese I wrote a few years back.

Notifications You must be signed in to change notification settings

ElioChen/Solution_for_Every_Gaussian_Error_Message

 
 

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

54 Commits
 
 
 
 

Repository files navigation

Solutions for Gaussian Error Messages (translation in progress ⏳)

This is an article containing reliable explanations and solutions for almost all error messages in the Gaussian quantum chemistry programs. This is a per-request English translation for a long blog in Chinese I wrote a few years back.

I wrote this as there is almost no reliable list of how to solve the errors. Similar online resources online are neither comprehensive nor accurate.

This blog is for people with basic knowledge of how Gaussian works, what's the basic structure of Gaussian files, and a basic understanding of quantum chemistry. So do not ask questions like "what do you mean by add XXXX keyword?".

Using (including part of) this blog in a commercial setting (or for any profitable action) is prohibited. Otherwise, feel free to use this material as long as you mention the source before the usage, and link to this article.

I'm not related to Gaussian Inc. in any way.

List of errors

General errors

Errors raised by specific links

Computer/OS related errors

CPU

File

Hard drive

Memory

Difficult cases

What does a Gaussian Bug look like

Abnormal situations that don't give an error message

Normal outputs that could be mistaken as errors

Other common questions

Explanation and solution for the errors

For ones who would like to read this in whole, the order of the following sections is different from the list of errors above. Use the anchor links if you would like to navigate to specific errors.

General errors

▶️ Severe Error Message # 2070 (Windows)

▶️ segmentation violation/segmentation fault (Linux)

➤ Example:

Linux

Error: segmentation violation
   rax 0000000000000000, rbx 000000000060c640, rcx ffffffffffffffff
   rdx 0000000000007f57, rsp 0000007fbffed948, rbp 0000007fbffed970
   rsi 000000000000000b, rdi 0000000000007f57, r8  0000002a9558af40
   r9  0000000000000000, r10 0000007fbffed801, r11 0000000000000206
   r12 0000000000615890, r13 0000000000640398, r14 0000000000640398
   r15 0000000000640398

Windows

➤ Explanation

This is just a message saying "Something is wrong" without any specific information.

➤ Solution

As there is no info on the specific error, any method you got online for "how to solve the Gaussian 2070 error" is nonsense. You should look at the (end of the) output file to see the specific reason for the error.

On the windows version of the Gaussian, the output might not update when the error happened, so you might see an empty window. In this case, you need to open the output file itself or wait some time to see the text in the main UI.

SCF Convergence

▶️ L502,Convergence failure

➤ Example:

 Cycle 128  Pass 1  IDiag  1:
 RMSU=  9.02D-07    CP:  1.00D+00  9.48D-01  3.00D+00  1.61D+00  1.30D+00
                    CP:  1.61D+00  1.46D+00
 E= -14914.8006510842     Delta-E=        0.000260457560 Rises=F Damp=F
 DIIS: error= 1.39D-02 at cycle 128 NSaved=  20.
 NSaved=20 IEnMin=16 EnMin= -14923.1431281454     IErMin= 8 ErrMin= 8.37D-03
 ErrMax= 1.39D-02  0.00D+00 EMaxC= 1.00D-01 BMatC= 8.92D-01 BMatP= 4.53D-01
 IDIUse=2 WtCom= 0.00D+00 WtEn= 1.00D+00
 EnCoef did    94 forward-backward iterations
 Coeff-En:   0.837D-01 0.000D+00 0.544D-01 0.287D+00 0.818D-03 0.000D+00
 Coeff-En:   0.000D+00 0.000D+00 0.149D+00 0.485D-05 0.228D-01 0.144D-03
 Coeff-En:   0.903D-02 0.129D-03 0.778D-04 0.393D+00 0.216D-03 0.313D-03
 Coeff-En:   0.237D-03 0.376D-03
 Coeff:      0.837D-01 0.000D+00 0.544D-01 0.287D+00 0.818D-03 0.000D+00
 Coeff:      0.000D+00 0.000D+00 0.149D+00 0.485D-05 0.228D-01 0.144D-03
 Coeff:      0.903D-02 0.129D-03 0.778D-04 0.393D+00 0.216D-03 0.313D-03
 Coeff:      0.237D-03 0.376D-03
 Gap=     0.015 Goal=   None    Shift=    0.000
 RMSDP=1.87D-06 MaxDP=3.06D-04 DE= 2.60D-04 OVMax= 8.02D-05

 >>>>>>>>>> Convergence criterion not met.
 Error on total polarization charges = ********
 SCF Done:  E(RPBE1PBE) =  -14914.8006511     A.U. after  129 cycles
            NFock=128  Conv=0.19D-05     -V/T= 4.2204
 KE= 4.631417055578D+03 PE=-4.911726716823D+04 EE= 7.213578297303D+03
 SMD-CDS (non-electrostatic) energy       (kcal/mol) =      -6.81
 (included in total energy above)
 Convergence failure -- run terminated.
 Error termination via Lnk1e in /gpfs/share/home/1501110295/g16/l502.exe at Tue Mar 31 23:16:51 2020.
 Job cpu time:       7 days 18 hours 47 minutes 16.6 seconds.
 Elapsed time:       0 days  5 hours 57 minutes 42.6 seconds.
 File lengths (MBytes):  RWF=   3744 Int=      0 D2E=      0 Chk=    176 Scr=     32

➤ Explanation

This is one of the most common problems in daily operations. The SCF (self-consistent field) iteration is not converged to below the threshold you designated (scf=conv=XX, default 1D-8) within the maximum number of cycles (scf=maxcyc=XX, default 128 cycles) you designated.

➤ Identify the problem

Before finding solutions, you need to identify specifically how is the convergence going. To do this, you need to replace the # to #p in your route card to show details of SCF iteration (I suggest always doing this, except for special jobs like ONIOM with a large molecule). You can plot the E = -XXX.XXXXX in each cycle, zoom in to the graph so that the changes in the later cycles are clearly visible (referred to as "the curve" later). From the curve, you may find a few scenarios, and you should deal with each one of them separately:

  1. [Fluctuation at initial stage]: The curve is fluctuating, and the convergence (given by NFock=128 Conv=XXXXD-XX in the output) is far from the threshold (usually larger than 10^-3 hartree). The fluctuation can be somewhat chaotic:

or can also be oscillating with a fixed period:

(please ignore the vertical axis, this is an output of my QM Monitor program, and the vertical axis does not linearly correspond to energy.)

  1. [Fluctuation at later stage]: same with above, just the convergence is closer to the threshold, like 10^-5 hartree or below

  2. [Monotonous decline]: The curve is monotonously declining (it's ok to have a few points as exceptions to the monotonous trend), but doesn't converge at maxcyc. This situation is exceedingly rare. I believe I've only encountered this less than 5 times in all my calculations, and interestingly, most of them involves lanthanide atoms.

  1. [SCF (curve shape irrelevant) with drastically different electronic energy for similar structures] (for example, during optimization). For example, the image and output below represent a job that ends with an SCF convergent failure. The last step has an SCF energy several hundred hartree higher than previous ones. Reviewing the optimization process (Image below, right), you can see the optimization kinda goes haywire (the forces and displacements don't look right) at around step 39. This is reflexed in the jump of his SCF energy (despite it has converged). Steps 50 and 56 (highlighted with arrow) also have energy tens of hartrees above the normal ones. This usually indicated a deeper problem than just SCF convergence problem and could be raised from dividing a very small number and results in a very large numerical error. The small number could be caused by the linear dependency of the basis sets or solvation cavity issues (see the section L502/L508,Inv3 failed in PCMMkU).

    Opt Step 37:  SCF Done:  E(RB-P86) =  -3825.13282335     A.U. after   12 cycles
    Opt Step 38:  SCF Done:  E(RB-P86) =  -3825.13719799     A.U. after   11 cycles
    --> Opt Step 39:  SCF Done:  E(RB-P86) =  -3825.10522704     A.U. after   13 cycles
    Opt Step 40:  SCF Done:  E(RB-P86) =  -3825.13602726     A.U. after   12 cycles
    Opt Step 41:  SCF Done:  E(RB-P86) =  -3825.13709759     A.U. after    8 cycles
    Opt Step 42:  SCF Done:  E(RB-P86) =  -3825.13718461     A.U. after    7 cycles
    Opt Step 43:  SCF Done:  E(RB-P86) =  -3825.13719706     A.U. after    1 cycles
    Opt Step 44:  SCF Done:  E(RB-P86) =  -3825.13719791     A.U. after    1 cycles
    Opt Step 45:  SCF Done:  E(RB-P86) =  -3825.13222362     A.U. after   12 cycles
    Opt Step 46:  SCF Done:  E(RB-P86) =  -3825.13968083     A.U. after   12 cycles
    Opt Step 47:  SCF Done:  E(RB-P86) =  -3825.13599635     A.U. after   13 cycles
    Opt Step 48:  SCF Done:  E(RB-P86) =  -3825.13662431     A.U. after   12 cycles
    Opt Step 49:  SCF Done:  E(RB-P86) =  -3825.13636208     A.U. after   12 cycles
    --> Opt Step 50:  SCF Done:  E(RB-P86) =  -3845.72449495     A.U. after   54 cycles
    Opt Step 51:  SCF Done:  E(RB-P86) =  -3825.13665330     A.U. after   20 cycles
    Opt Step 52:  SCF Done:  E(RB-P86) =  -3825.13647203     A.U. after   23 cycles
    Opt Step 53:  SCF Done:  E(RB-P86) =  -3825.13594059     A.U. after   26 cycles
    Opt Step 54:  SCF Done:  E(RB-P86) =  -3825.13482223     A.U. after   26 cycles
    Opt Step 55:  SCF Done:  E(RB-P86) =  -3825.13408631     A.U. after   35 cycles
    --> Opt Step 56:  SCF Done:  E(RB-P86) =  -3888.91182114     A.U. after   73 cycles
    Opt Step 57:  SCF Done:  E(RB-P86) =  -3825.13665328     A.U. after   22 cycles
    Opt Step 58:  SCF Done:  E(RB-P86) =  -3825.13647358     A.U. after   25 cycles
    Opt Step 59:  SCF Done:  E(RB-P86) =  -3825.13595614     A.U. after   26 cycles
    Opt Step 60:  SCF Done:  E(RB-P86) =  -3825.13490799     A.U. after   27 cycles
    Opt Step 61:  SCF Done:  E(RB-P86) =  -3825.13382018     A.U. after   30 cycles
    --> Opt Step 62:  SCF Done:  E(RB-P86) =  -4207.17296836     A.U. after  129 cycles
    

➤ Solution

First, there is almost no keyword that can systematically increase the chance of convergence in everyday calculation. I suggest against adding any additional SCF-related keywords to your input file (i.e., keep it as the default).

(Many of the following content in this section referenced this Chinese post by Sobereva.)

  1. For [Monotonous decline] ONLY, use scf=maxcyc=XX (choose a number larger than 128) to increase the maximum number of cycles. As I mentioned before, this is a rare situation, and this keyword is heavily abused on the internet. Use this keyword in any other situations will only cause you unnecessary CPU times. And it is a sign of ignorance if you see someone write this option as his "default option" (i.e. write this to all the input files).

  2. For [Fluctuation at initial stage], first try one or a combination of several of the keywords below:

    • SCF=NoVarAcc: In L502 output, you can see this sentence Initial convergence to 1.0D-05 achieved. Increase integral accuracy. (example below), which means before this step, as the SCF is far from convergent, Gaussian will use a lower integration accuracy to save time. However, this may cause SCF fluctuation, in that case, disabling it may solve the problem. Adding this keyword will almost never cause a convergent SCF to become divergent, but it will likely cause an increase in the computation time.

           Coeff:      0.113D+00 0.365D+00 0.534D+00
           Gap=     0.125 Goal=   None    Shift=    0.000
           RMSDP=9.29D-06 MaxDP=3.04D-03 DE=-3.18D-03 OVMax= 3.90D-03
      
           Cycle  10  Pass 0  IDiag  1:
           RMSU=  6.05D-06    CP:  9.87D-01  6.50D-01  1.65D-01  6.63D-02  6.41D-01
                              CP:  6.58D-01  7.67D-01  8.17D-01  7.24D-01
           E= -13297.5461801855     Delta-E=       -0.000664594365 Rises=F Damp=F
           DIIS: error= 2.05D-04 at cycle  10 NSaved=  10.
           NSaved=10 IEnMin=10 EnMin= -13297.5461801855     IErMin=10 ErrMin= 2.05D-04
           ErrMax= 2.05D-04  0.00D+00 EMaxC= 1.00D-01 BMatC= 1.37D-04 BMatP= 8.51D-04
           IDIUse=3 WtCom= 9.98D-01 WtEn= 2.05D-03
           Coeff-Com: -0.513D-03 0.181D-02-0.342D-02 0.137D-02-0.404D-02-0.120D-01
           Coeff-Com:  0.518D-02 0.669D-01 0.176D+00 0.769D+00
           Coeff-En:   0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00
           Coeff-En:   0.000D+00 0.000D+00 0.000D+00 0.100D+01
           Coeff:     -0.511D-03 0.180D-02-0.341D-02 0.137D-02-0.404D-02-0.120D-01
           Coeff:      0.517D-02 0.668D-01 0.176D+00 0.769D+00
           Gap=     0.123 Goal=   None    Shift=    0.000
           RMSDP=4.95D-06 MaxDP=6.20D-04 DE=-6.65D-04 OVMax= 1.54D-03
      
           --> Initial convergence to 1.0D-05 achieved.  Increase integral accuracy.
           Cycle  11  Pass 1  IDiag  1:
           E= -13297.5469692442     Delta-E=       -0.000789058700 Rises=F Damp=F
           DIIS: error= 2.66D-04 at cycle   1 NSaved=   1.
           NSaved= 1 IEnMin= 1 EnMin= -13297.5469692442     IErMin= 1 ErrMin= 2.66D-04
           ErrMax= 2.66D-04  0.00D+00 EMaxC= 1.00D-01 BMatC= 3.14D-04 BMatP= 3.14D-04
           IDIUse=3 WtCom= 9.97D-01 WtEn= 2.66D-03
           Coeff-Com:  0.100D+01
           Coeff-En:   0.100D+01
           Coeff:      0.100D+01
           Gap=     0.117 Goal=   None    Shift=    0.000
           RMSDP=4.95D-06 MaxDP=6.20D-04 DE=-7.89D-04 OVMax= 4.47D-03
      
           Cycle  12  Pass 1  IDiag  1:
           RMSU=  7.87D-05    CP:  1.00D+00
           E= -13297.5470667379     Delta-E=       -0.000097493670 Rises=F Damp=F
      
    • SCF=NoIncFock: Incremental Fock matrix build is an acceleration technique where the Fock matrix is computed recursively using the difference density of the last 2 SCF cycles. This could drastically lower the scaling of the Fock matrix build. However, this may cause fluctuation, especially when diffuse functions are used. To solve this, use SCF=NoIncFock turns off incremental Fock matrix build. This will result in higher cost of each SCF cycle.

    • IOp(5/37=N): By default, Gaussian will rebuild the Fock matrix in full for every 20 cycles of SCF iteration: IOp(5/37)=20 by default. This is shown as Fock matrices will be formed incrementally for 20 cycles. and Restarting incremental Fock formation. in the Gaussian output (example of output below). You can select an N<20, and involve the IOp(5/37=N) option to form the Fock matrix every N SCF cycles. The relationship between SCF=NoIncFock and IOp(5/37=N) are similar to opt=CalcAll and opt=Recalc=N.

       Requested convergence on             energy=1.00D-06.
       No special actions if energy rises.
       --> Fock matrices will be formed incrementally for  20 cycles.
      
       Cycle   1  Pass 1  IDiag  1:
       FoFJK:  IHMeth= 1 ICntrl=       0 DoSepK=F KAlg= 0 I1Cent=           0 FoldK=F
       IRaf= 990000000 NMat=       1 IRICut=       1 DoRegI=T DoRafI=F ISym2E= 0 IDoP0=0 IntGTp=1.
       FoFCou: FMM=T IPFlag=           0 FMFlag=      100000 FMFlg1=        2001
       ...
       ...
       ...
       Coeff:     -0.634D-05 0.481D-03-0.104D-02 0.243D-02-0.167D-02 0.429D-02
       Coeff:      0.217D-02 0.976D-02-0.371D-01-0.324D-01 0.394D-02 0.132D+00
       Coeff:      0.254D+00 0.663D+00
       Gap=     0.029 Goal=   None    Shift=    0.000
       RMSDP=1.13D-03 MaxDP=7.30D-02 DE=-4.18D-04 OVMax= 1.76D-02
      
       Cycle  21  Pass 1  IDiag  1:
       --> Restarting incremental Fock formation.
       E= -1315.53915604307     Delta-E=       -0.000105907567 Rises=F Damp=F
       DIIS: error= 1.06D-04 at cycle  21 NSaved=  20.
       NSaved=20 IEnMin=20 EnMin= -1315.53915604307     IErMin=20 ErrMin= 1.06D-04
       ...
       ...
       ...
       Coeff:      0.627D-01 0.170D+00-0.268D+00-0.291D+00-0.212D+01 0.175D+01
       Coeff:     -0.254D+01 0.241D+01-0.587D+00 0.111D+01 0.166D+01-0.269D+01
       Coeff:     -0.561D+01 0.530D+01-0.209D+01 0.473D+01
       Gap=     0.052 Goal=   None    Shift=    0.000
       RMSDP=7.12D-04 MaxDP=8.59D-02 DE=-5.43D-03 OVMax= 4.24D-03
      
       Cycle  41  Pass 1  IDiag  1:
       --> Restarting incremental Fock formation.
       E= -5974.52121524936     Delta-E=       -0.000886211828 Rises=F Damp=F
       DIIS: error= 8.43D-03 at cycle  41 NSaved=  17.
       NSaved=17 IEnMin=11 EnMin= -5976.91129255226     IErMin=10 ErrMin= 6.51D-03
       ```
      
    • SCF=vshift=N:

      • Use level shift to virtually increase the HOMO-LUMO gap. SCF=vshift=N shifts orbital energies by N*0.001 hartree. You can usually use a value of 200~500. This does NOT affect the final converged results. It's most effective to systems with small gaps, like (multicore) transition metal complexes. Use this if other method fails.
    • SCF=Fermi:

      • A black box mixed method uses Fermi broadening, damping and level shifting dynamically. Use this if other method fails.
  • Make sure you have a qualitatively correct wavefunction:

    • There are situations where the guess you provided to the program is qualitatively incorrect.

    • For example, the guess for the system of [HSO3-]...[NO2·] could have an incorrect fragment charge population and corresponds to [HSO3·]...[NO2-], causing SCF problems (even if you got a converged wavefunction, the result would be meaningless for your purpose). Similar situation could arise from spin population: the most common problem is a closed shell guess was given for a spin polarized system; more complicated situation could be the spin population for dual-transition metal complexes.

    • For these situations, the most effective way to solve this is a fragmented guess, where a two-step Gaussian job is called, and the first step only generates a guess (where the wavefunction for each fragment is converged, then put together without further optimization).

    • The following input file is for the calculation of [HSO3-]...[NO2·]. In the second step, guess=read cause Gaussian to read the wavefunction from the chk file generated above; guess=always let Gaussian read the guess.chk at every optimization step, this can sometimes prevent the drift-away from the wavefunction you want during optimization. pop=always asks Gaussian to calculate Mulliken population for every step of the optimization, so that you can check that the wavefunction is always correct.

      %rwf=/path/to/your/temp/temp.rwf
      %nosave
      %chk=/path/to/your/work_dir/guess.chk
      #p PBE1PBE def2SV guess=fragment=2
      
      Generate fragmented wavefunction
      
      -1 2 -1 1 0 2                               <-- The charge/multiplicity for the whole system, for [HSO3-] (Fragment 1), and for [NO2·] (Fragment 2) respectively
      H(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      S(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      N(fragment=2)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=2)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=2)  XX.XXXX  XX.XXXX  XX.XXXX
      
      --link1--
      %rwf=/path/to/your/temp/temp.rwf
      %nosave
      %oldchk=/path/to/your/work_dir/guess.chk
      %chk=/path/to/your/work_dir/calculation.chk
      #p PBE1PBE def2SV guess=(read,always) opt freq pop=always
      
      Read fragmented wavefunction for the calculation
      
      -1 2                                        <-- The charge/multiplicity for the whole system only
      H(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      S(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=1)  XX.XXXX  XX.XXXX  XX.XXXX
      N(fragment=2)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=2)  XX.XXXX  XX.XXXX  XX.XXXX
      O(fragment=2)  XX.XXXX  XX.XXXX  XX.XXXX
      
      
    • Be aware there are more complicated situations, but it is out of scope of this document.

  • Use a converged wavefunction as initial guess:

    • A converged wavefunction (if selected wisely) could be a much better guess than the initial guess generated by Harris or other guess method, avoiding fluctuation at initial stage. Possible choice includes:

    • If you are doing calculation with a large basis set (especially when there are diffuse functions), use converged function at smaller basis set (or without diffuse functions). You can use a looser SCF convergent threshold in the initial guess step. The following file generates a guess with def2-SV(P) basis set, and use it as initial guess for calculation at ma-TZVP level:

      %rwf=/path/to/your/temp/temp.rwf
      %nosave
      %chk=/path/to/your/work_dir/guess.chk
      #p PBE1PBE def2SV SCF=conv=5
      
      Generate guess wavefunction
      
      0 1
      [coordinates]
      
      --link1--
      %rwf=/path/to/your/temp/temp.rwf
      %nosave
      %oldchk=/path/to/your/work_dir/guess.chk
      %chk=/path/to/your/work_dir/calculation.chk
      #p PBE1PBE gen guess=(read) opt freq
      
      Read converged wavefunction at smaller basis set for the calculation
      
      0 1
      [coordinates]
      
      @/path/to/your/ma-TZVP.basis
      
    • If you are doing a calculation with an anion (-1 1), use the cation (2-electron less, 1 1) as the guess. This could work with neutral molecules, where you use the (2 1) state as the guess, but that has a lesser effect than the anion to cation case.

    • If you are doing calculation with a radical (0 2), use the closed shell cation as the guess (1 1).

    • If you know a similar but different geometry that can results in a converged SCF, use that as a guess. This usually works when a bond is being broken, while a longer or shorter bond distance could work. Also be aware of potential polarized singlet problem. And sometimes just randomly change the geometry also works (if that is the case, be vigilant whether there are problems about solvent cavity)

    • If you are running a calculation with solvation, use vacuum SCF wavefunction as a guess.

    • If you are running a calculation with pure functional, or a hybrid functional with low Hartree-Fock percentage, try to converge a wavefunction with Hartree-Fock, and use it as the guess. HF is easier to converge as it has a larger (overestimate) gap.

    • If you are calculating an ROHF (ROKS) state, use UHF (UKS) wavefunction as a guess.

  • Use a different guess method:

    • Use keywords guess=huckel or guess=INDO. However, I found this to be rarely useful.
  1. For [Fluctuation at later stage]:
    • Int=UltraFine:

      • 【When it is useful】For methods that use a grid-based numerical integration (DFT only), insufficient grid quality could cause small numerical error that's just larger than the convergent threshold of SCF (and for the same reason, this is not likely to help fluctuations at early stage of SCF). The Minnesota functionals are kinda notorious for their higher requirement for grid quality. Thus, this is the prime recommendation if that is the case. I have examples where this keyword help with other functionals, but it's rarer than the Minnesota functionals.
      • 【G16】Note that since Gaussian 16, the default option has already become Int=Ultrafine. So unless you specifically lowered it to Int=Fine, there is no need to write Int=Ultrafine explicitly.
      • 【Comparablity】Energies at different integration grid is not comparable. If you raise the integration grid for one structure, you need to re-run all the calculations for which you are getting a difference in energy with that structure. One exception is that you can optimize the structure with Int=UltraFine grid, but do an Int=Find grid single point calculation to get the correct electronic energy, and add the later electronic energy with the thermodynamic correction at Int=UltraFine grid (This is because theintegrationn grid is not likely to greatly affect the PES). You need to clearly state what integration grid you used in the Supporting Information of your publication (for repeatability).
      • 【About SuperFine】There is a higher grid Int=SuperFine, however Int=UltraFine is usually good enough, and it is unlikely that SuperFine grid will solve SCF divergence, unless you are requiring some uncommon SCF threshold.
    • SCF=NoIncFock

    • IOp(5/37=N)

  2. Danger zone
    • SCF=QC, SCF=XQC, SCF=YQC:
      • Asks the use of quadratically convergent method (for SCF=XQC, the QC is invoked only if ordinary SCF failed after (by default) 32 steps). SCF=QC has a significant better chance of convergence. However, it's (1) far more expensive, (2) more likely to result in a local minimum or unconverged wavefunction (do sanity check with wavefunctional analysis, or at least do a stable calculation, especially with ONLY SCF=QC could converge but not any other method).
    • SCF=conver=N:
      • Use coarser criteria for SCF convergence. Default is 8, so choose a number smaller than 8. The specific meaning is to converge the RMS of the density matrix to below 10-N and converge the maximum change of the density matrix and the change of electronic energy to below 10-N+2. For single point calculations, you can set it to scf=conver=6 with reasonable result. However, do not use this in optimization or frequency missions. As if you converge the SCF wavefunction to 10-6, the gradient would be only be converged to at about 10-3 ~ 10-4. The default optimization convergence criteria for gradient is on the order of 10-4, which means you may never be able to converge a geometry optimization.
    • IOp(5/13=1):
      • This IOp suppress the checking for SCF convergence. The program will carry out the following reaction no matter what's the situation of the SCF convergence. This should almost NEVER be used. The result is likely to be qualitatively wrong. If anyone suggests this as a solution to a L502 error, or use this option in his "default input file template", he should not be trusted to give any advice for Gaussian, and all of this computational result should be checked before interpretation.
      • Some people say you can use it in a geometry optimization, as if one of the steps is not converged, the next step could, and finally leads to a reasonable result. However, this is nonsense. If the SCF is not converged, the gradient has a much larger error than the energy, which means the following optimization step is basically a random move. If this solves the convergence problem, you may as well do this by yourself. Also, one may well be forgot to check EVERY optimization mission that the SCF is indeed converged in the last step, and report a meaningless result.

▶️ L508,Convergence failure ⏳

Density matrix breaks symmetry, PCut= 6.91D-03
Density matrix has no symmetry -- integrals replicated.
Iteration  80 EE= -1377.03506721338     Delta-E=       -0.001137227868 Grad=4.268D-02
Gradient too large for Newton-Raphson -- use scaled steepest descent instead.
Convergence failure.
Error termination via Lnk1e in C:\G09W\l508.exe at Mon Jun 26 19:36:04 2017.

成因:用 QC 方法时SCF不收敛

解决方法:去掉SCF=qc、SCF=xqc、SCF=yqc 关键词进行计算。如果继而出现 L502, Convergence Failure问题,参照上文解决,只不过绕开涉及 scf=qc 的方案即可。

Post-HF Convergence

▶️ L913,*MAX. CYCLES* ⏳

CCSD、CCSD(T):

Iteration Nr.   50
**********************
DD1Dir will call FoFMem   1 times, MxPair=         2
NAB=     1 NAA=     0 NBB=     0.
Norm of the A-vectors is  4.0812478D-03 conv= 1.00D-05.
RLE energy=       -0.0380878830
DE(Corr)= -0.37790236E-01 E(CORR)=     -1.1531999202     Delta=-1.55D-03
NORM(A)=   0.10071894D+01
*************
*MAX. CYCLES*
*************
Largest amplitude= 4.58D-02
Error termination via Lnk1e in C:\G09W\l913.exe at Sun Mar 05 10:12:08 2017.
Job cpu time:       0 days  0 hours  0 minutes  1.0 seconds.
File lengths (MBytes):  RWF=     22 Int=      0 D2E=      0 Chk=      1 Scr=      1

CISD

Iteration Nr.  50
**********************
DD1Dir will call FoFMem   1 times, MxPair=         2
NAB=     1 NAA=     0 NBB=     0.
Norm of the A-vectors is  4.3269814D+01 conv= 1.00D-05.
RLE energy=       -9.6293060152
DE(CI)=       5.1143514     E(CI)=        4.3978787385   
NORM(A)=   0.36740722D+02
SIZE-CONSISTENCY CORRECTION:
S.C.C.=    0.00000000D+00        E(CI,SIZE)= -0.71597715082D+00
*************
*MAX. CYCLES*
*************
***************************************************************
Dominant configurations:
***********************
Spin Case        I    J    A    B          Value
    AA            1         3           -0.354624D+01
......
   ABAB           1    1   23   23       0.147097D+00
Largest amplitude= 1.84D+01
Error termination via Lnk1e in C:\G09W\l913.exe at Sat Jun 11 03:43:10 2016.
Job cpu time:       0 days  0 hours  0 minutes  2.0 seconds.
File lengths (MBytes):  RWF=     23 Int=      0 D2E=      0 Chk=      1 Scr=      1

成因:高斯用某种迭代方法求解 CISD 和 CCSD(CCSD(T)),而该迭代过程也会不收敛,在默认循环(50步)或设定的循环内未收敛即会报此错误。(注意CISD报出 Max Cycles 错误之后还会吐一堆组态,真正报错在它上面)

解决: 如下">>>"标出的一行所示,计算过程中会输出每步迭代后的CCSD / CISD相关能(DE(corr)),CCSD还会输出本轮迭代的能量变化(Delta,即本轮的 DE(corr) 减去上一轮的 DE(corr)),需观察前者是否有稳定的收敛至某一值、后者是否逐渐趋近于0,以判断是否有收敛趋势。

**********************
DD1Dir will call FoFMem   1 times, MxPair=       380
NAB=   190 NAA=     0 NBB=     0.
Norm of the A-vectors is  1.2424967D-05 conv= 1.00D-05.
RLE energy=       -1.1041652032
>>> DE(Corr)=  -1.1041652     E(CORR)=     -236.55647322     Delta=-7.56D-08
NORM(A)=   0.11814199D+01
Iteration Nr.  12
**********************

如果上述能量有收敛趋势,增加CCSD、CCSD(T)、CISD迭代循环的最大值,默认为50,应大于此值,语法如 CCSD(T) 改为 CCSD(T, maxcyc=100) 若能量在震荡,应首先检查结构、参考态是否合理,如超过稳定点的单重态双自由基是否用了对称性破缺的初猜;另外可尝试微调结构、换基组 如果最终震荡的幅度很小,例如能量的默认收敛限为1E-7,最终Delta一直在 nE-7的水平,到不了E-8,可以用诸如 CCSD(T,conver=6) 将收敛阈值提高至 1E-6。(注意如果是在热力学组合方法中遇到这个问题,就不应该写CCSD(T,conver=6),而应该只写CCSD(conver=6)。

Geometry Optimization Convergence

▶️ L9999,Optimization stopped. ⏳

Error termination via Lnk1e in l9999.exe

此报错需查看输出文件中的额外信息,自结尾处向上检索“Optimization stop”可见如下几种情况

Optimization stopped.
-- Wrong number of Negative eigenvalues: Desired= 1 Actual= 4
-- Flag reset to prevent archiving.

成因:Gaussian会在过渡态优化时自动检查虚频的数量,若不为1即终止。这一检查没有必要 解决:将NoEigenTest关键词加入到opt的选项中,如 opt=(TS, CalcFC, NoEigenTest)

Optimization stopped.
-- Number of steps exceeded, NStep= 100
-- Flag reset to prevent archiving.

成因:几何优化未能在指定步数内收敛 解决:打开GaussView,File-Open,按照如下设置: image

然后点击Result-Optimization,可以看到优化时的能量变化曲线。观察优化曲线应纵向放大,方法为用鼠标反复选择纵轴方向仅包含最末几个点、但横轴方向包含所有点的“宽而矮”的矩形范围,直到能看清最末几个点的变化情况为止。根据图形判断是否震荡。如下图,需纵向放大至右边才能看清确实震荡了。(注意图中偶尔有突然的、无“周期性的”某1~2个点的spike是正常现象,不是震荡,下方小幅度但有“周期性”的才是震荡。)

左侧看不出是否震荡。需纵向放大至右侧才能看出。 image image

如果没有震荡:opt=maxcycle=当前步数的2~3倍。可以读取当前最末/最优的一个构象做为新的初猜。注意检查构象、轨迹是否合理,如果结构已经不合理了,则应适当调整初猜后再做(尤其是找过渡态的时候)。

如果震荡了,看:http://sobereva.com/164。

文中的众多解决方法中,我个人倾向于找到能量较低、4个收敛标准离收敛限较近,且结构合理的某点为新的初猜,并首先尝试在 opt 选项中加上(MaxStep=5, NoTrustUpdate, GDIIS) 选项;如果有计算频率的资源,可再加上opt=calcfc选项。特别仅对于对称性较高的结构,如T、C3 等群的复杂分子,若反复调整关键词后仍然震荡,可以尝试构建或破坏具有对称性的初猜结构,分别结合opt=Cartesian 优化,如 opt=(MaxStep=5,NoTrustUpdate,Cartesian), 有时有效果。 对于某些结构,SMD溶剂模型可能带来数值问题,此时震荡不是几何优化算法造成的,不论如何调整opt的关键字也不会解决。此时可以在气相、或SMD的其他溶剂空腔选项下优化(若需与其他结构级别统一,可以在所得结构的基础上用与其他结构相同的溶剂模型算个单点)。 注意上一段仅是我个人首先尝试的建议,不要仅仅试了这一个不行就又把问题拿出来问。应该继续尝试 http://sobereva.com/164 的其他办法。

▶️ L123,Max corrector steps exceded ⏳

▶️ L123,GS2 Optimization Failure. ⏳

Delta-x Convergence NOT Met
    Maximum number of corrector steps exceded.
    Error termination via Lnk1e in l123.exe at Sun Sep  6 07:30:18 2015.
GS2 Optimization Failure.
    Error termination via Lnk1e in l123.exe

成因:Gaussian的默认IRC算法HPC为了使曲线平滑需做重校正,但常有校正步不收敛的问题。换用GS2算法时,其限制性优化有时也不收敛。

解决:增加IRC=LQA选项,例如IRC=(calcfc,LQA)。

Geometry Manipulation

▶️ L103,Error in internal coordinate system ⏳

▶️ L103,Linear angle in Bend; Linear angle in Tors ⏳

▶️ L103,FormBX had a problem ⏳

Bend failed for angle     1 -    11 -     3
    Tors failed for dihedral     9 -     1 -    11 -     3
    Tors failed for dihedral    10 -     1 -    11 -     3
    Tors failed for dihedral    12 -     1 -    11 -     3
    Tors failed for dihedral    14 -     3 -    11 -     1
    Tors failed for dihedral    17 -     3 -    11 -     1
    FormBX had a problem.
    Error termination via Lnk1e in l103.exe
(Enter /home/gauuser/g09/l103.exe)

GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
Berny optimization.
Using GEDIIS/GDIIS optimizer.
Linear angle in Bend.
Error termination via Lnk1e in /home/gauuser/g09/l103.exe at Sun Feb  5 05:25:49 2017.
NTrRot=    -1 NTRed=   798 NAtoms=    66 NSkip=   606 IsLin=F
Error in internal coordinate system.
Error termination via Lnk1e in l103.exe
Berny optimization.
Using GEDIIS/GDIIS optimizer.
Linear angle in Tors.
Error termination via Lnk1e in /home/program/g09/l103.exe at Thu Jun 25 22:29:51 2015.

成因:内坐标有其内在的限制,遇到优化过程中遇到有几个原子正好排成直线等情况可能出现此问题 解决: (0) 检查轨迹、当前结构、冻结(柔性扫描)变量设置是否合理 (1) opt=cartesian 这一方法可从原理上彻底解决此问题,但 opt=cartesian 在多数情况下会增加优化至相应极小点所需的步数,如果体系不非常耗时可以直接用其至优化收敛。 如果体系接近自己计算能力的极限,可以用 opt=cartesian 计算两三步之后终止计算,保存坐标后重新写一个输入文件,换回默认的opt方法。 另外注意在使用modredundant冻结、扫描了坐标的时候,不能与opt=cartesian合用 (2) 有时直接保存最末结构后,重新opt也能解决此问题,Gaussian 实际上会自动为接近直线的原子添加一些 Linear Bend,但并不总是有效。 (3) 部分情况下,人为增加一 Linear bend 冗余内坐标,输入为在冗余内坐标段落加入如 L 1 2 3 -1 B 的表述,并启用Modredundant,其中 1 2 3 为 “Bend failed for angle” 的三个原子,如以上图报错来说,若原输入文件为:

#p B3LYP/genecp opt freq

Title

0 1
[原子坐标们]

C O H
TZVP
****
Co
SDD
****

Co
SDD

则新的输入文件应为(注意共3处不同):

#p B3LYP/genecp opt=modredundant freq

Title

0 1
[自出错的输出文件中提取的、出错前一步结构的原子坐标们]

L 1 11 3 -1 B

C O H
TZVP
****
Co
SDD
****

Co
SDD

其中 “[自出错的输出文件中提取的、出错前一步结构的原子坐标们]” 可以用GaussView打开之前出错的输出文件,随便另存为一个gjf输入文件,用文本编辑器打开,把坐标部分拷出来。 此处用 genecp 并添加相应自定义基组、赝势的段落是为了说明输入文件中各段落的顺序。 如果你不知道我在说什么,可阅读:http://sobereva.com/60 完整的输入文件段落顺序可见于:http://sobereva.com/g09/m_input.htm

▶️ L103,New curvilinear step not converged. Error imposing constraints. ⏳

Iteration 96 RMS(Cart)=  0.00000206 RMS(Int)=  0.00542712
Iteration 97 RMS(Cart)=  0.00000193 RMS(Int)=  0.00542766
Iteration 98 RMS(Cart)=  0.00000180 RMS(Int)=  0.00542817
Iteration 99 RMS(Cart)=  0.00000169 RMS(Int)=  0.00542865
Iteration100 RMS(Cart)=  0.00000158 RMS(Int)=  0.00542909
New curvilinear step not converged.
Error imposing constraints
Error termination via Lnk1e in /home/gauuser/g09/l103.exe at Mon Jul 17 10:21:34 2017.

成因:进行需要做限制性优化的任务(如opt=modredundant 的 F 和 S、QST2等),Optimizer 不知道在当前限制条件下该把结构初猜摆成什么样子。

解决: 如果用了QST2: 换用TS(Berny)方法找过渡态,QST2不是什么好方法

如果在柔性扫描中遇到问题: (1)如果是在从上一个扫描点(如键长为1.7A)、至下一个扫描点的过程中产生的问题(如键长为1.8A),则可以手动将初猜结构中的相应键长调整到下一点的键长(1.8A),从这一点开始后续扫描。最后把结果拼接起来就可以了 (2)扫描过程中频繁出现此问题有可能是“步子扯得太大”,可以把步长调小 (3)如果有很快速、合适的方法进行预扫描(如PM6-D3)可以先在相应级别下做扫描,产生每一点的结构,然后用多个限制性优化(F)代替扫描(S),此时由于并行化是100%的,只要预优化合理、效率还比柔性扫描高。

如果在限制性优化中(或柔性扫描的某一点中间)遇到问题: (1)稍微调整结构重新来 (2)偶尔将结构重新保存为输入文件、直接接着运行问题也可以解决。原因不明。可以尝试

Nonsense Input File (Formatting / Structure / Route)

▶️ L101,End of file in ZSymb. ⏳

▶️ L301,End of file reading basis center. ⏳

▶️ L301,EOF while reading ECP pointer card. ⏳

   C -1.21995   2.13345   0.
   End of file in ZSymb.
   Error termination via Lnk1e in l101.exe
End of file reading basis center.
Error termination via Lnk1e in /home/gauuser/g09/l301.exe at Wed Mar  1 20:50:57 2017.
======================================================================================================
                                       Pseudopotential Parameters
======================================================================================================
  Center     Atomic      Valence      Angular      Power
  Number     Number     Electrons     Momentum     of R      Exponent        Coefficient   SO-Coeffient
======================================================================================================
EOF while reading ECP pointer card.
Error termination via Lnk1e in /home/gauuser/g09/l301.exe at Wed Mar  1 20:43:41 2017.

成因:Gaussian需要读到换行表示的段落结尾才认为相应输入段落结束,若当前段落最后一个字符后,没有空行,Gaussian会读入一个EOF(End of File)特殊字符,不明白啥意思就跪了。这两个问题的提示不同是因为输入文件中分别以原子坐标段落(没有gen)、混合基组定义(gen)、和赝势定义(genecp)的段落为末尾,从而在读取不同信息时报错了。

解决:输入文件最后多打几个空行(应总是这么做,没坏处)。

▶️ L101,WANTED A XXX AS INPUT. FOUND AN XXX AS INPUT. ⏳

Wanted an integer as input.
Found a string as input.
H                0. 0. 0.                                                      
?
Error termination via Lnk1e in C:\G09W\l101.exe at Sun Mar 26 17:48:52 2017.
Job cpu time:       0 days  0 hours  0 minutes  0.0 seconds.
File lengths (MBytes):  RWF=      5 Int=      0 D2E=      0 Chk=      1 Scr=      1

成因:输入文件格式错误。高斯会检查特定段落的变量类型是否正确,诸如电荷、多重度应为整数(如 “1”),坐标应为浮点数(如 “1.0”),还有文本类型诸如Title段落等等,若不正确即报此错误。

解决:根据提示具体分析错误在哪儿。如故意制造上面报错使用的输入文件为:

#p B3LYP/6-31G(d)

Title

H 0. 0. 0.

由报错提示“H 0. 0. 0. ”可知Gaussian在这一行想读到其他东西,而不是坐标,仔细思考就知道是缺少了电荷和多重度,都是整型变量,读到了“H”是文本型变量,故而报错,修改正确即可:

#p B3LYP/6-31G(d)

Title

0 2
H 0. 0. 0.

▶️ L101,Input Error Input Error Input Error Input Error Input Error Input Error ⏳

-------------------
Title Card Required
-------------------
Symbolic Z-matrix:
Charge =  0 Multiplicity = 1


Input Error Input Error Input Error Input Error Input Error Input Error

There are no atoms in this input structure !

Please fix the molecule specification section of your input and try again.

Input Error Input Error Input Error Input Error Input Error Input Error

Error termination via Lnk1e in C:\G09W\l101.exe at Wed Jan 16 20:26:00 2019.
Job cpu time:       0 days  0 hours  0 minutes  0.0 seconds.
File lengths (MBytes):  RWF=      5 Int=      0 D2E=      0 Chk=      1 Scr=      1

成因:Gaussian没有读到结构声明

解决:检查:(1)是不是真没写结构;(2)是不是忘了写geom=check, geom=allcheck;(3)是不是多了/少了空行,使高斯读入了空段落。

▶️ L1,Illegal IType or MSType generated by parse. ⏳

----------
#p sp freq
----------
Illegal IType or MSType generated by parse.
Error termination via Lnk1e in C:\G09W\l1.exe at Thu Dec 07 13:58:19 2017.
Job cpu time:       0 days  0 hours  0 minutes  0.0 seconds.
File lengths (MBytes):  RWF=      1 Int=      0 D2E=      0 Chk=      1 Scr=      1

成因:撰写输入文件时缺乏基本常识,将不搭界的、相互矛盾的计算任务写在同一个计算路径里。例如上面的输入告诉高斯要算单点、同时又说要算频率,高斯懵逼了。 解决:别犯傻。(不要自己乱组合, 不要让 sp 和 freq,force 和 freq 等矛盾的任务设定在一条计算路径里出现)

▶️ L1,QPErr --- A syntax error was detected in the input line. ⏳

------------------
#p M06-2X/6-31G(d)
------------------
QPErr --- A syntax error was detected in the input line.
#p M06-2X/6-31G(d)
       '
Last state= "GCL"
TCursr=      3656 LCursr=         7
Error termination via Lnk1e in C:\G09W\l1.exe at Wed Mar 06 14:19:36 2019.
Job cpu time:       0 days  0 hours  0 minutes  0.0 seconds.

成因:关键词或语法错误,错误的位置由下一行的 “ ’ ” 标记给出了 解决:按着手册写关键词,纠正关键词或语法错误。例如上面的例子中,程序标记出M06-2X关键字有问题,看手册得知高斯中正确的泛函写法应为M062X,修改后即可。

▶️ L301,The combination of multiplicity X and XXX electrons is impossible. ⏳

(Enter C:\G09W\l301.exe)
Standard basis: 6-31G(d) (6D, 7F)
Ernie: Thresh=  0.10000D-02 Tol=  0.10000D-05 Strict=F.
The combination of multiplicity 1 and     1 electrons is impossible.
Error termination via Lnk1e in C:\G09W\l301.exe at Sun Mar 26 18:00:32 2017.
Job cpu time:       0 days  0 hours  0 minutes  1.0 seconds.
File lengths (MBytes):  RWF=      5 Int=      0 D2E=      0 Chk=      1 Scr=      1

成因:完全没正确设置/理解自旋多重度或电荷。自旋多重度定义如下,即 Alpha电子数 - Beta电子数 + 1,一个电子要么是 alpha 电子,要么是 beta 电子。故可知偶数电子的体系,自旋多重度只能是奇数(1,3,5,7...),反之亦然。 解决:正确设置自旋多重度、电荷、结构。 应检查诸如:自由基的自旋多重度习惯性的写了1;正负离子体系忘了改分子净电荷数;建模时结构中多了、少了氢原子之类的问题。 另外需注意,GaussView虽然看似会“自动判断”自旋多重度,但只是提供最低可能的自旋多重度而已(即偶数电子给1,奇数电子给2),应该用结构化学的知识自己判断体系可能的自旋多重度可能是多少(如闭壳层分子、单重态双自由基为1,一般有机自由基为2,金属配合物应该用配体场理论判断等等),如果实在不确定(如对某些过渡金属配合物),可以在不同的自旋多重度下优化分子,看哪种自旋多重度的能量最低。

▶️ L202,Problem with the distance matrix. ⏳

Small interatomic distances encountered:      6     1     7     2     8     3     9     4    10     5
Problem with the distance matrix.
Error termination via Lnk1e in C:\G09W\l202.exe at Wed Jan 16 20:33:11 2019.

成因:有至少两个原子间距离太近,很有可能距离是0。例如上述的报错中的一串数字写明了是:6号原子和1号原子之间太近,7号和2号太近,8号和3号太近,etc 解决:找到相应原子调整分子结构,或者删去多余的原子。一种出现这种情况的原因是在GaussView中复制粘贴了一遍自己的结构而不自知,所以在GaussView里看着特别正常,但打开原子列表就可以发现问题了。

▶️ L202,Atoms too close. ⏳

Small interatomic distances encountered:
     2    1 5.00D-02
Atoms too close.
Error termination via Lnk1e in C:\G09W\l202.exe at Wed Jan 16 20:38:12 2019.

成因:原子间距离太近,小于高斯认为正常结构应有的阈值(但不很接近于0)。例如上述例子中报错是2号原子和1号原子之间距离太近,2号和1号的距离是 0.05A。 解决:一般情况下,应该检查结构,确认其正常,拉远相应的原子;如果是在进行扫描键长的任务、扫描曲线某一端刻意要求某两个原子离得很近,可以使用Geom=NoCrowd关闭这个检查,但注意即使是在扫描里,这么近的地方一般也没什么意义;在其他任务里用 Geom=NoCrowd 就更掩耳盗铃了。如果是在例如优化等过程中走了几步变成这样了,一定是级别(方法、基组)或其他设置有低级错误,诸如某个原子未给定基组之类,应仔细检查。

▶️ L301,Atomic number out of range for XXX basis set. ⏳

Rotational constants (GHZ):      0.0817250      0.0474806      0.0408748
Standard basis: 6-31G(d) (6D, 7F)
Atomic number out of range for 6-31G basis set.
Error termination via Lnk1e in /data2/G09/g09/l301.exe at Fri Mar 24 20:58:27 2017.

成因:使用内置基组进行计算时,分子中含有基组不支持(超出基组定义范围)的元素。如上面的案例中,计算使用了6-31G(d)基组,其只支持H~Kr;但结构中含有Pt原子,故报错为Pt超出基组定义范围。 解决:选择合适的基组(及赝势),按照规定格式输入,可参考:

谈谈量子化学中基组的选择 http://sobereva.com/336 详解Gaussian中混合基组、自定义基组和赝势基组的输入 http://sobereva.com/60爬虫并处理后的离线EMSL基组库 http://bbs.keinsci.com/forum.php ... hlight=%BB%F9%D7%E9 在线基组和赝势数据库一览 http://sobereva.com/309

▶️ L301,R6DS8: Unable to choose the S8 parameter ⏳

IExCor=  408 DFT=T Ex=B Corr=PW91 ExCW=0 ScaHFX=  0.000000
ScaDFX=  1.000000  1.000000  1.000000  1.000000 ScalE2=  1.000000  1.000000
IRadAn=      0 IRanWt=     -1 IRanGd=            0 ICorTp=0 IEmpDi=141
NAtoms=    1 NActive=    1 NUniq=    1 SFac= 1.00D+00 NAtFMM=   60 NAOKFM=F Big=F
Integral buffers will be    262144 words long.
Raffenetti 2 integral format.
Two-electron integral symmetry is turned on.
R6DS8: Unable to choose the S8 parameter, IExCor=  408 IXCFnc=  0 ScaHFX=  0.000000 IDFTD=4
Error termination via Lnk1e in C:\G09W\l301.exe at Wed Jul 12 20:20:24 2017.

成因:使用 empiricaldispersion=GD3 或 empiricaldispersion=GD3BJ 时程序未存有相应泛函的D3色散校正参数。 解决:(1) 确认在D3或D3BJ中选择的正确,例如米尼苏达系列泛函没有D3-BJ参数,不能用empiricaldispersion=GD3BJ,只能用empiricaldispersion=GD3 (2) 换一个有色散校正的泛函 (3) 自定义色散校正参数(往往来源于较新的文献),可以参考《Gaussian中非内置的理论方法和泛函的用法》(http://sobereva.com/344)一文中2.2节的做法,以及手册中的相应部分 (4) 不用色散校正

▶️ L301,R6DRCv: No RCov radius available for IA=XX ⏳

▶️ L301,R6DC6: No C6 coefficient available for IA=XX ⏳

Integral buffers will be    262144 words long.
Raffenetti 2 integral format.
Two-electron integral symmetry is turned on.
R6DC6: No C6 coefficient available for IA= 96
Error termination via Lnk1e in C:\G09W\l301.exe at Sat Mar 17 13:56:16 2018.
Integral buffers will be    262144 words long.
Raffenetti 2 integral format.
Two-electron integral symmetry is turned on.
R6DRCv: No RCov radius available for IA=                  96
Error termination via Lnk1e in C:\G09W\l301.exe at Sat Mar 17 13:57:31 2018.

成因:使用了色散校正(empiricaldispersion=GD3,empiricaldispersion=GD3BJ,empiricaldispersion=GD2,以及使用了自带DFT-D的泛函比如wB97xD),且程序中没有相应元素的参数。上文的两个例子分别是在使用DFT-D2和DFT-D3时,分子中有锔原子(IA=XX 即为相应的元素序数)时的报错。例如对 DFT-D3 来说,原文只为 1-94 号元素(H-Pu)做了参数化,对 96 号元素使用 DFT-D3 当然会缺参数。

解决:不用色散校正。

▶️ L301,No solvent atoms in DisRep. ⏳

Solvent              : n,n-DiMethylFormamide, Eps=  37.219000 Eps(inf)=   2.046330
                        RSolv=   0.000000 Ang.
------------------------------------------------------------------------------
Warning! Inconsistent VMol and RSolv for this solvent, using VMol=    0.00 Ang**3.
No solvent atoms in DisRep.
Error termination via Lnk1e in /home/igors/g09/l301.exe at Fri May 19 02:16:27 2017

成因:

示例出错的输入文件如下(下面的乙酸溶剂没有缺乏相应非极性参数):

# opt() freq m062x/6-311G(2df,p) SCRF(CPCM, solvent=AceticAcid, read)

Title Card Required

0 1
O                  0.00000000    0.00000000    0.11813800
H                  0.00000000    0.75681000   -0.47255200
H                  0.00000000   -0.75681000   -0.47255200

dis
rep
cav

即使用PCM溶剂模型时使用 “Dis Rep Cav” 考虑了非极性部分,而所使用的相应溶剂在高斯中并没有定义非极性部分的参数。高斯列表支持的溶剂当中(G09 或 G16 有约180种)只有 G03 年代就支持的约20种有相应参数支持以这样的方式考虑非极性部分(可参考 G03 的说明书),其他均会报此错误。

解决:高斯或许可能有方法自定义 Dis Rep Cav 的参数,但未在说明书中给出。可以用SMD模型考虑非极性部分(对自定义溶剂可见:http://bbs.keinsci.com/forum.php?mod=viewthread&tid=6683);或是不考虑非极性部分。

▶️ L602,GetVDW: no radius for atom XX atomic number XX. ⏳

(Enter /home/software/g09D01/g09/l602.exe)
FitSet:  NAtFit=    76 NAtPot=    76 NAtFrz=     0 MDM=    80 TotChg=   0.00000
Merz-Kollman atomic radii used.
GetVDW:  no radius for atom   2 atomic number  29.
Error termination via Lnk1e in /home/software/g09D01/g09/l602.exe at Mon Dec 26 13:57:24 2016.
Job cpu time:       0 days  4 hours  6 minutes 15.9 seconds.
File lengths (MBytes):  RWF=    856 Int=      0 D2E=      0 Chk=     31 Scr=      1

原因:使用pop=CHELPG 拟合静电势时没有内置相应元素的半径(上面例子中是第29号元素Cu)。 解决:pop里用readradii,输入文件末尾写上元素名和指定的半径(一般用范德华半径,可以查得),例如此例对铜来说应写如下的关键词:

#p (其他关键词) pop=(CHELPG,ReadRadii)

title

坐标们

C H N P Cl 0
6-31G*
****
Cu 0
lanl2dz
****

Cu 0
lanl2dz

Cu 1.4

这里用赝势是为了体现各段落的顺序关系

▶️ L1002,No func 3rd derivs with XXX. ⏳

Isotropic polarizability for W=    0.000000       12.56 Bohr**3.
No func 3rd derivs with HSE.
Error termination via Lnk1e in D:\Program Files (x86)\g09\G09W\l1002.exe at Sun Mar 26 11:50:16 2017.
Job cpu time:       0 days  0 hours  0 minutes  5.0 seconds.
File lengths (MBytes):  RWF=      5 Int=      0 D2E=      0 Chk=      1 Scr=      1

成因:此任务要求Gaussian用LC-wPBE泛函计算超极化率,polar任务计算超极化率需要三阶解析导数且默认认为DFT方法都有三阶解析导数,但可见于Gaussian手册对DFT方法的说明,少数泛函没有:

AVAILABILITY

Energies, analytic gradients, and analytic frequencies; ADMP calculations.

Third order properties such as hyperpolarizabilities and Raman intensities are not available for functionals for which third derivatives are not implemented: the exchange functionals Gill96, P (Perdew86), BRx, PKZB, TPSS, wPBEh and PBEh; the correlation functionals PKZB and TPSS; the hybrid functionals OHSE1PBE and OHSE2PBE.

解决:使用具有三阶导数的方法;或令程序用数值方法计算三阶解析导数 具体至使用LC-wPBE计算极化率的问题: (1)用 polar=Numerical 计算数值方法的极化率α (2)可以使用默认的Polar关键词,虽然其在计算超极化率β时会因为LC-wPBE没有三阶导数而报错,但报错前极化率α实际上已经输出来了(如下面段落“>>>”所示的行)。

FullF1:  Do perturbations    1 to     3.
>>> SCF Polarizability for W=    0.000000:
>>> 1             2             3
>>> 1  0.311786D+02
>>> 2  0.214915D+01  0.273009D+02
>>> 3 -0.264577D-04  0.351431D-04  0.193900D+02
Isotropic polarizability for W=    0.000000       25.96 Bohr**3.
No func 3rd derivs with HSE.
Error termination via Lnk1e in D:\Program Files (x86)\G09\G09W\l1002.exe at Thu Apr 13 08:57:29 2017.
Job cpu time:       0 days  0 hours  0 minutes 14.0 seconds.

▶️ L103,NRF ne Abs NRFX ⏳

Cartesian Forces: Max 0.101934913 RMS 0. 010471640
Leave Link 716 at Sun Jul 3 00:08:55 2016, MaxMem= 1073741824 cpu:0.7
(Enter /home/export/base/fanyu/apps/g09/1103.exe)

GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
Berny optimization.
Internal Forces: Max 0.077268082 RMS 0.018267038
Search for a local minimum.
Step number 1 out of a maximum of 202
All quantities printed in internal units (Hartrees-Bohrs-Radians)
NRF ne Abs NRFX
Error termination via Lnkle in /home/export/base/fanyu/apps/g09/1103.exe at Sun Jul 3 00:08:57 2016.

成因:在ONIOM计算中冻结不合理。此例中其冻结了所有 Low layer 而放开了所有 High layer。 对 Link atom 进行处理时,由于Link atom的种类在两层中不同(如 Si-O 键打断,用 Si-H 封端),两层中的键长应当不同;所以 Link atom 和它在 High layer 里连接的原子要么都冻结、要么都放开。

解决:

  1. link atom和其连接的原子用相同方式处理
  2. 对于中等大小的体系(<300个原子),可以使用Opt=NoMicro来取消MicroIteration

详细讨论见:http://bbs.keinsci.com/forum.php?mod=viewthread&tid=3671

Computer / OS Related Problems

▶️ Error: illegal instruction, illegal opcode ⏳

这一错误不在输出文件中报错,而是在 Linux 终端中显示,常见样式如下:

Error: illegal instruction, illegal opcode
   rax 000000001b23ef84, rbx 0000000000730438, rcx 000000001b23ef84
   rdx 00007fffa89b0b10, rsp 00007fffa89b0b40, rbp 00007fffa89b0b50
   rsi 00000000007c29f0, rdi 00007fffa89b0b10, r8  0000003e2d78fee8
   r9  0000000000000001, r10 00007fffa89b0880, r11 0000000000000202
   r12 00007fffa89c7878, r13 00007fffa89c7848, r14 0000000000000000
   r15 0000000000000000
  /lib64/libpthread.so.0() [0x3e2e00f710]
  /home/gauss/g09/l1.exe() [0x4a72e1]
  /home/gauss/g09/l1.exe() [0x41d295]
  /home/gauss/g09/l1.exe() [0x403745]

输出文件中没有有效信息,仅输出初始命令

Entering Gaussian System, Link 0=g09
Input=CO.com
Output=CO.log
Initial command:
/home/gauss/g09/l1.exe "/home/gauss/g09/tmp/Gau-4412.inp" -scrdir="/home/gauss/g09/tmp/"

成因: Gaussian的可执行文件有不同的指令集版本,不同指令集的CPU需购买不同版本的Gaussian程序(具体可见:http://gaussian.com/g16/g16_plat.pdf)

自己的CPU支持什么指令集可以在Linux上使用命令“cat /proc/cpuinfo”查看,若Gaussian的版本是“AVX2-enabled”,则仅当 “/proc/cpuinfo” 文件的 flags 字段中支持AVX2时才能使用,否则应该更换其他指令集版本的Gaussian程序,如 Legacy 即 pre-SSE4.2 版本

解决: 更换Gaussian的可执行文件的指令集版本,如将Gaussian 09 E.01 X86_64(AVX2-enabled) (兼容性差、但性能高) 更换为 Gaussian 09 E.01 X86_64(“Legacy”-pre-SSE4.2)(下兼容性好,但性能稍差) 更换机器、使用最新的CPU

▶️ Internal input file was deleted! ⏳

Leave Link  703 at Wed Mar  1 20:23:53 2017, MaxMem= 18983419904 cpu:       101.8
Internal input file was deleted!
Error termination via Lnk1e at Wed Mar  1 20:23:53 2017.

该报错总发生于某一Link刚刚结束、新Link尚未开始之时,但具体发生在哪一Link并不固定。

成因:Gaussian似乎需要一个完整的(没有“@”)、当前步的输入文件来完成Link之间的交接,这个临时文件存储于 $GAUSS_SCRDIR 环境变量下,文件名和路径都会在输出文件中写明,如下面的任务中,该临时文件就是“/home/gauuser/g09/scratch/Gau-16535.inp”。这个报错的原因是这个文件somehow在运行之中被删掉了。

Initial command:
/home/gauuser/g09/l1.exe "/home/gauuser/g09/scratch/Gau-16535.inp" -scrdir="/home/gauuser/g09/scratch/"
Entering Link 1 = /home/gauuser/g09/l1.exe PID=     16536.

解决:重算(或视任务类型,从当前步骤Restart),并尝试指认这个临时文件被删除的原因并解决它。遇到过的情况包括:使用了一些奇怪的脚本、如在某任务 A 结束后清除临时文件夹,但却误将任务 B 的临时文件也清除了;队列控制、文件系统等问题;当然还包括自己、或使用同一机器的其他人手贱→_→。

▶️ open-new-file ⏳

Gaussian 09:  IA32W-G09RevD.01 24-Apr-2013
                05-Nov-2017
******************************************
%chk=C:\gtest\C6H5NO2.chk
fname=C:\gtest\C6H5NO2.chk fd = -1

open-new-file

原因:高斯无法创建、访问 %chk(或其他Link0指令)指向的文件(本例中是 C:\gtest\C6H5NO2.chk) 解决: 确认相应目录确实存在(此例中是 C:\gtest 目录)。 确认高斯对其有读写权限(Windows下赋予高斯(G09W.exe)管理员权限,Linux 下对调用高斯的用户调整相应文件夹所有权与权限)。 确认相应文件没有正被占用。

▶️ L1,ntrex1 ⏳

%mem=130000MB
%chk=D:\Gaussian\ORCaussian_test\Parallel\ORCA_derv\calcfc_opt_TS2\TS.chk
ntrex1

成因、解决:Gaussian无法访问指定的路径,例如在Linux上写了Windows下的路径,改正确了即可。

▶️ Error Message # 2066. Can't create file Temp input file 'gxx.inp' ⏳

image

成因:Gaussian无法创建临时文件

解决:确认Gaussian 09W File-Preference 中的 Scratch Path 目录存在;确认Gaussian对该目录有读写权限(可以尝试修改至系统盘之外的其他磁盘;或给予 Gaussian09w.exe 管理员权限);确认磁盘未满。

image

▶️ Erroneous write ⏳

Erroneous write. Write -1 instead of 4096.
    fd = 4
    orig len = 4096 left = 4096
    g_write

成因:硬盘不够 解决: image

▶️ g_write ⏳

▶️ L914,XXXXXX words are not enough for AIAXAO. ⏳

▶️ L1002,XXXXXX words are not enough for AIAXAO. ⏳

▶️ needs more words of memory ⏳

▶️ Not enough memory to run at all ⏳

▶️ Out-of-memory error in routine XXX ⏳

   Out-of-memory error in routine FoFDir-SEAll (IEnd=             1 MxCore=   -1974194264)
    Use %mem=1915MW to provide the minimum amount of memory required to complete this step.
    Error termination via Lnk1e in l914.exe at Tue Mar 08 15:32:05 2016.
Iteration     1 Dimension    60 NMult    60
Cannot handle 2e integral symmetry, ISym2E=1.
CISAX needs   2533811 more words of memory.
Error termination via Lnk1e in l914.exe
Not enough memory to run at all:  LenERI=     2254414
Increase memory by     1247333 words.
Error termination via Lnk1e in d:\gaussian03\l804.exe
GetIJB would need an additional   555279372 words of memory to use all  12 processors.
    JobTyp=2 Pass  1:  I=  41 to  74 NPSUse=  3 ParTrn=F ParDer=T DoDerP=T.
Generate precomputed XC quadrature information.
          Solving linear equations simultaneously, MaxMat=      72.
AlAXAO:  NMat=   690 NPMax=    1 NPMax1=    0 MaxMat=    72.
    33484188 words are not enough for AlAXAO.
Error termination via Lnk1e in /cluster/gaussian-09.C.01-amd64-sse4a//g09/l1002.exe at Tue Jan  6 15:17:39 2015.
Convergence on wavefunction:    0.001000000000000
Iteration     1 Dimension   128 NMult     0 NNew    128
AlAXAO:  NMat=   128 NPMax=    1 NPMax1=    0 MaxMat=     0 Max3X=F.
    33548910 words are not enough for AlAXAO.
Error termination via Lnk1e in /public/home/qsh/g09/l914.exe at Tue Jul 11 13:23:34 2017.

成因:提示都非常直接,内存不够 解决:用%mem设置增加内存即可。 对某些情况,如“need an additional words of memory to use all XX processors.”,若没有足够的内存提供,可以减少运行核数(%nprocshared)。若仍然不够,加内存条 or 降低计算机别

▶️ galloc: could not allocate memory. ⏳

Leave Link    1 at Fri Mar  8 01:49:58 2013, MaxMem=17716740096 cpu:       0.8
galloc:  could not allocate memory.

成因:输入文件中准许Gaussian程序使用的内存超过了机器当前剩余的内存量,故Gaussian无法想系统申请到想要数量的内存而报错。 解决:将输入文件中 %mem 命令后的内存量设置为当前系统剩余的内存量,并留下足够的余量(防止运行途中剩余内存减少而又出现这个问题)。在Windows下查看剩余内存可用任务管理器;在 Linux 下用 top 命令查看;若同时运行多个任务,也需额外注意内存分配问题。

Others

▶️ L502/L1002,Inaccurate quadrature in CalDSu. ⏳

Inaccurate quadrature in CalDSu.
Error termination via Lnk1e in l502.exe
Inaccurate quadrature in CalDSu.
Error termination via Lnk1e in l1002.exe

首先应检查: (1)基组/赝势是否合理,是否有例如赝势基组没写赝势之类的低级错误 (2)结构(或是轨迹)是否合理

如果均没有问题,则按下面的PPT解决:

# #### ▶️ L114,ERROR IN INITNF. NUMBER OF VARIABLES ( 0) INCORRECT (SHOULD BE BETWEEN 1 AND 50) ⏳
NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-NEF-
NUMERICAL EIGENVECTOR FOLLOWING MINIMUM SEARCH
INITIALIZATION PASS


************************************************
** ERROR IN INITNF. NUMBER OF VARIABLES (  0) **
**   INCORRECT (SHOULD BE BETWEEN 1 AND 50)   **
************************************************


Error termination via Lnk1e in C:\G09W\l114.exe at Sun Mar 26 16:57:35 2017.
Job cpu time:       0 days  0 hours  0 minutes  0.0 seconds.
File lengths (MBytes):  RWF=      5 Int=      0 D2E=      0 Chk=      1 Scr=      1

成因:使用无解析梯度的方法做优化而为明确声明变量,如CCSD(T)方法 解决: 更换方法,除非极特殊的目的,否则用CCSD(T)之类的方法优化总是没有必要的,如果体系很小、计算量很充裕,用CCSD优化已经足够了。 按下面 Sobereva 的 PPT 处理: image image image

▶️ L801,Excessive mixing of core and valence orbitals. ⏳

Largest valence mixing into a core orbital is  5.06D-01
Largest core mixing into a valence orbital is  4.89D-01
Excessive mixing of frozen core and valence orbitals.
Error termination via Lnk1e in l801.exe

成因:Gaussian在post-HF、双杂化、TD等运算时会做默认的冻芯,不将其考虑到激发行列式中,这有时会导致问题 解决: (1)最简单的解决方法是在理论方法的关键词中加入(Full)选项,如做MP2计算出问题时对应用MP2(Full)代替、做CCSD(T)时用CCSD(T,full)代替,使程序不使用冻芯近似。 (2)做其他的冻芯设定、可见Gaussian手册对Frozen Core的说明(http://sobereva.com/g09/k_fc.htm) (3)IOp(8/11=1)可防止在检查得到冻芯不合理时即终止任务,而仅将这个警告打印出来、接着进行后续运算。但切记自己检查结果的合理性。

▶️ L801, Fatal Problem: The smallest alpha delta epsilon is XXXXXX ⏳

(Enter /home/bjwang/g09/l801.exe)
Range of M.O.s used for correlation:     1   264
NBasis=   264 NAE=   114 NBE=   114 NFC=     0 NFV=     0
NROrb=    264 NOA=   114 NOB=   114 NVA=   150 NVB=   150

**** Fatal Problem: The smallest alpha delta epsilon is -0.78537370D-01

Error termination via Lnk1e in /home/bjwang/g09/l801.exe at Mon Jan 31 18:17:13 2011.

成因:在TDDFT、post-HF 中做积分变换时涉及分子轨道的能量差为分母的项,若分子轨道Gap太小(“alpha delta epsilon” 即两个α轨道的△ε),分母接近0,得到一个大数会导致数值不稳定性问题。具体参看本文后面关于“Warning!!: The smallest alpha delta epsilon is...”的章节。

解决:由于未能找到可重复的算例,以下解决方案仅为推测、未经完全证实,仅供尝试 (1)IOp(8/11=1)可防止在出现此警告时即终止任务,仅将这个警告打印出来、接着进行后续运算。但切记自己检查结果的合理性。 (2)换泛函,更换HF成分较高的泛函可能能解决这个问题(增大Gap)

▶️ L502/L508,Inv3 failed in PCMMkU ⏳

(Enter /home/g16/l502.exe)
Integral symmetry usage will be decided dynamically.
Closed shell SCF:
Using DIIS extrapolation, IDIIS=  1040.
NGot= 18983288832 LenX= 18982597916 LenY= 18982263254
Requested convergence on RMS density matrix=1.00D-08 within 128 cycles.
Requested convergence on MAX density matrix=1.00D-06.
Requested convergence on             energy=1.00D-06.
No special actions if energy rises.
Fock matrices will be formed incrementally for  20 cycles.

Cycle   1  Pass 1  IDiag  1:
FoFJK:  IHMeth= 1 ICntrl=       0 DoSepK=F KAlg= 0 I1Cent=           0 FoldK=F
IRaf= 990000000 NMat=       1 IRICut=       1 DoRegI=T DoRafI=F ISym2E= 0 IDoP0=0 IntGTp=1.
FoFCou: FMM=T IPFlag=           0 FMFlag=      100000 FMFlg1=        2001
         NFxFlg=           0 DoJE=F BraDBF=F KetDBF=F FulRan=T
         wScrn=  0.000000 ICntrl=         0 IOpCl=  0 I1Cent=           0 NGrid=           0
         NMat0=    1 NMatS0=      1 NMatT0=    0 NMatD0=    1 NMtDS0=    0 NMtDT0=    0
Symmetry not used in FoFCou.
FMM levels:  10  Number of levels for PrismC:   9
Inv3:  Mode=1 IEnd=    55186563.
Iteration    1 A*A^-1 deviation from unit magnitude is 2.18D-11 for   1454.
Iteration    1 A*A^-1 deviation from orthogonality  is 1.36D-11 for   4245   1454.
Iteration    1 A^-1*A deviation from unit magnitude is 2.55D-11 for   1454.
Iteration    1 A^-1*A deviation from orthogonality  is 2.25D-10 for   4227   3067.
Iteration    2 A*A^-1 deviation from unit magnitude is 7.03D-07 for   4227.
Iteration    2 A*A^-1 deviation from orthogonality  is 8.89D-07 for   2879   1454.
Iteration    2 A^-1*A deviation from unit magnitude is 1.00D-06 for   1454.
Iteration    2 A^-1*A deviation from orthogonality  is 9.53D-07 for   4227   1454.
...
...
Iteration   10 A*A^-1 deviation from unit magnitude is 7.97D-07 for   2943.
Iteration   10 A*A^-1 deviation from orthogonality  is 1.08D-06 for   2943   1454.
Iteration   10 A^-1*A deviation from unit magnitude is 3.30D-07 for   4227.
Iteration   10 A^-1*A deviation from orthogonality  is 6.08D-07 for   4227   2061.
Inv3 failed in PCMMkU.
Error termination via Lnk1e in /home/g16/l502.exe at Mon Sep  4 11:33:20 2017.
Job cpu time:       0 days  0 hours 18 minutes 30.5 seconds.
Elapsed time:       0 days  0 hours  0 minutes 42.9 seconds.

这一报错发生于SCF的第一次迭代(含正常的 L502 和 QC 的 L508)。

Gaussian Help 回复大意如下(译自英文):

这个问题与溶剂模型的溶质空腔有关。高斯中为构建SMD溶质空腔的默认选项是为中等大小的溶质分子设置的(注:默认是范德华表面 surface=vdW),它是由以原子为中心、一定的原子半径的一些球构成的。其中用于构建空腔的默认原子半径较小,对较大的溶质分子会造成溶剂空腔长得像“奶酪”,有很多数学形式上暴露于溶剂中的孔洞、沟槽,但在实际体系中这些位置并不能接触溶剂。这会造成上述为 PCM 矩阵求逆迭代过程中不收敛的问题,也是你另一封邮件中提到SCF能量异常的本质原因(我同时发送了两封邮件,下文将提到)。

理想情况下,使用溶剂排斥表面(Solvent-excluding surface, SES)是更好的选择。SES 首先以原子位置为中心,构建出范德华表面;之后再用一些原点不在原子上的球面进行平滑,就消除了上述不能被溶剂够到的孔洞、沟槽,解决此问题。但其劣势在于这个方法做出的表面不随构象的改变而连续变化,造成势能面不连续、不平滑,难以进行几何优化。

另一个可能是用溶剂可及表面(Solvent-accessible surface,SAS)构建空腔。溶剂可及表面和 vdW 表面一样都是由以原子为中心的球构成的,但其半径是[原子半径+溶剂半径](注,这个时候溶剂半径才有用,默认选项下自定义溶剂时设半径没什么用)。SAS 表面没有“奶酪问题”,且势能面连续。但其问题是半径太大,空腔扣的太多,相比 vdW 和 SES 表面会低估溶剂化效应。

一个较好的解决方法是用 SAS 表面做优化,然后用 SES 表面做单点。或者在 SAS 表面做优化后再用其结果为初猜做 SES 表面下的单点(以期待 SES 在极小点附近可能没有不连续性)。这两种做法都比在气相优化的结果好一些。

image

关键词:

更换溶剂表面的类型用 scrf=read 关键词启用溶剂化控制的段落,并在相应段落声明 surface=xxx,如 surface=SAS。 示例如下:

#p pbe1pbe/genecp opt freq
empiricaldispersion=gd3bj
scrf=(smd,solvent=water,read)

title

0 1
C                   1.629576   -0.586781    1.988383
....
H                   3.604870   -0.932068   -2.049232

[Basis set definition]
...

[Pseudopotential definition]
...

surface=sas

对 surface=SES,还需添加 AddSph 选项。即将上述文件中的 “surface=sas” 换为 “surface=SES AddSph”

如果需要,可以用 geom=allcheck 等语句接续下一步计算切换为其他表面做单点/优化。

这个问题也有可能不报错,但造成 SCF 能量明显异常。这个错误比较隐蔽,当出现能量明显不对的时候应该检查。

例如我在一次优化中发现有两次 SCF 迭代后的能量(用>>>标出)低了好几十个 Hartree,明显是不对的,SCF 迭代所需的次数也明显较多。这也是上面所述的溶剂空腔问题。

        Line 68580:  SCF Done:  E(RB-P86) =  -3825.13599635     A.U. after   13 cycles
        Line 69843:  SCF Done:  E(RB-P86) =  -3825.13662431     A.U. after   12 cycles
        Line 71105:  SCF Done:  E(RB-P86) =  -3825.13636208     A.U. after   12 cycles
>>> Line 73188:  SCF Done:  E(RB-P86) =  -3845.72449495     A.U. after   54 cycles
        Line 74770:  SCF Done:  E(RB-P86) =  -3825.13665330     A.U. after   20 cycles
        Line 76261:  SCF Done:  E(RB-P86) =  -3825.13647203     A.U. after   23 cycles
        Line 77804:  SCF Done:  E(RB-P86) =  -3825.13594059     A.U. after   26 cycles
        Line 79343:  SCF Done:  E(RB-P86) =  -3825.13482223     A.U. after   26 cycles
        Line 81098:  SCF Done:  E(RB-P86) =  -3825.13408631     A.U. after   35 cycles
>>> Line 83672:  SCF Done:  E(RB-P86) =  -3888.91182114     A.U. after   73 cycles
        Line 85289:  SCF Done:  E(RB-P86) =  -3825.13665328     A.U. after   22 cycles
        Line 86819:  SCF Done:  E(RB-P86) =  -3825.13647358     A.U. after   25 cycles
        Line 88367:  SCF Done:  E(RB-P86) =  -3825.13595614     A.U. after   26 cycles
        Line 89942:  SCF Done:  E(RB-P86) =  -3825.13490799     A.U. after   27 cycles

▶️ L123,GetHes: LRWHes > LHess! ⏳

    ******** Start new reaction path calculation ********
RCFC Option Requested - Data Read From Chk File:
"/home/gauuser/Gaussian/PavB_Rad_Cyc/TS3/To_IM3__To_IM3_3[Complete_PBE1PBE]_Step
6[IRC_Forward].chk"
   Energy From Chk =  -2089.1954828

Current Structure is TS -> form Hessian eigenvectors.
                            Diagonalizing Hessian.
WARNING: NO IMAGINARY FREQUENCIES AT TS!
Supplied step size of   0.1000 bohr.
    Integration on MW PES will use step size of   0.2173 sqrt(amu)*bohr.
Point Number:   0          Path Number:   1


...


Leave Link  716 at Wed Jun 28 18:41:00 2017, MaxMem= 18666618880 cpu:         0.5
(Enter /home/gauuser/g09/l123.exe)
IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC
          ******** Start new reaction path calculation ********
LRWHes =   0  LHess  = ***
GetHes: LRWHes > LHess!
Error termination via Lnk1e in /home/gauuser/g09/l123.exe at Wed Jun 28 18:41:00 2017.
Job cpu time:       0 days  0 hours 56 minutes 39.2 seconds.

成因:此错误较为罕见,未能获得更多的错误样本。仅有的一个例子中,经分析是由于 1.读取了前一步 opt=TS freq 的Hessian,2. 该Hessian无虚频(注意前面有“WARNING: NO IMAGINARY FREQUENCIES AT TS!”提示)。 解决:检查过渡态的频率是否正常、只有一个虚频。若虚频数不对、首先解决过渡态优化的问题再说;如果没有该问题可上传文件讨论。

Difficult cases

▶️ NtrErr Called from XXXXXX. ⏳

Structure from the checkpoint file:  "sideroidide_CS_1_opt.chk"
FileIO operation on non-existent file.
FileIO: IOper= 2 IFilNo(1)=  -584 Len=          36 IPos=           0 Q=   47012767336112


dumping /fiocom/, unit = 1 NFiles =    33 SizExt =    524288 WInBlk =       512
                   defal = T LstWrd =     1206272 FType=2 FMxFil=10000

Number              0           501           502           503           507           511
Base           216576         43520         72192        180736        181760        196608
End           1206272         44520         76297        180821        181925        197835
End1          1206272         44544         76800        181248        182272        198144
Wr Pntr        216576         43520         72192        180736        181760        196608
Rd Pntr        216576         44520         72192        180736        181760        196608
Length         989696          1000          4105            85           165          1227

Number            551           552           561           562           575           579
Base           213504        212480        214016        205824        183808        212992
End            213542        212501        214017        212016        196499        213016
End1           214016        212992        214528        212480        196608        213504
Wr Pntr        213504        212480        214016        205824        183808        212992
Rd Pntr        213504        212480        214016        205824        183808        212992
Length             38            21             1          6192         12691            24

Number            598           603           665           672           674           698
Base            76800        216064        181248        183296        182784        215040
End             76802        216065        181684        183506        182924        215112
End1            77312        216576        181760        183808        183296        215552
Wr Pntr         76800        216064        181248        183296        182784        215040
Rd Pntr         76800        216064        181248        183296        182784        215040
Length              2             1           436           210           140            72

Number            700           701           730           761           801           989
Base           215552        198144        182272        214528        180224         44544
End            215842        205553        182514        214529        180230         64544
End1           216064        205824        182784        215040        180736         65024
Wr Pntr        215842        198144        182272        214528        180224         44544
Rd Pntr        215552        198144        182502        214528        180224         44544
Length            290          7409           242             1             6         20000

Number            991           992           993           994           995           996
Base            65536         65024         43008         40448         42496         41472
End             72098         65033         43208         40478         42516         41672
End1            72192         65536         43520         40960         43008         41984
Wr Pntr         65536         65024         43008         40448         42496         41472
Rd Pntr         72098         65033         43208         40478         42516         41672
Length           6562             9           200            30            20           200

Number            997           998           999
Base            41984         40960         77312
End             42274         41160        179816
End1            42496         41472        180224
Wr Pntr         41984         40960         77312
Rd Pntr         42274         41160         79816
Length            290           200        102504


dumping /fiocom/, unit = 2 NFiles =     1 SizExt =         0 WInBlk =       512
                   defal = F LstWrd =       65536 FType=2 FMxFil=10000

Number              0
Base            40448
End             65536
End1            65536
Wr Pntr         40448
Rd Pntr         40448
Length          25088


dumping /fiocom/, unit = 3 NFiles =     1 SizExt =    524288 WInBlk =       512
                   defal = T LstWrd =       65536 FType=2 FMxFil=10000

Number              0
Base            40448
End             65536
End1            65536
Wr Pntr         40448
Rd Pntr         40448
Length          25088
FileIO operation on non-existent file.
FileIO: IOper= 2 IFilNo(1)=  -584 Len=          36 IPos=           0 Q=   47012767336112
Error termination in NtrErr:
NtrErr Called from FileIO.
%chk=temp.chk
Bad file opened by FileIO:  Unit=2 I=  2 FPrev=40960 FCur=    0.
FileIO: IOper= 9 IFilNo(1)=     2 Len=           0 IPos=           0 Q=        135419884


dumping /fiocom/, unit = 1 NFiles =     1 SizExt =    524288 WInBlk =       512
                   defal = T LstWrd =       65536 FType=2 FMxFil=10000

Number           0
Base         20480
End          65536
End1         65536
Wr Pntr      20480
Rd Pntr      20480
Length       45056
Error termination in NtrErr:
NtrErr Called from FileIO.

NtrErr Called from XXXXXX的报错成因非常复杂,且真实报错离末尾较远,初学者常只贴 Error termination in NtrErr: NtrErr Called from XXXXXX. 乃至最后一屏幕等,实际上在其上的一堆数字之上,例如最开始的案例的实际报错是,后面一堆乱码一样的文字对用户没啥用:

Structure from the checkpoint file:  "sideroidide_CS_1_opt.chk"
FileIO operation on non-existent file.

此时就很有可能是chk文件出了什么问题。应首先检查chk文件是否存在、可读,例如可以检查: (1)路径是否正确、完整(如路径存在中文、存在空格、存在括号、没写扩展名等低级错误) (2)使用相对路径时,“相对”的不对,如两次运行程序时CWD(PWD)不同,导致不能找到同一个文件,可以将完整路径(Windows从盘符开始,Linux从/开始) (3)文件权限问题。在Linux上检查路径存在、相应文件(文件夹)有正确的读写权限。在Windows上,不要将临时文件存储在需要特殊权限的目录下,或者给Gaussian程序管理员权限后再运行。 (4)文件是否存在,例如我见过某些脚本在计算完成后会自动将chk文件删除、或将chk文件转化为fchk文件,此时当然会出错。 (5)如果实在找不到本质、应当先去除所有需要用到chk文件做输入的因素后重试,例如放弃使用geom=allcheck, guess=read之类的关键词从头开始算试试,可解决一大部分报错。这在本文一开始的部分也提到过。

注意类似因chk文件问题产生的报错有时并不直接,如下面这个案例中:

******************************************
Gaussian 09:  ES64L-G09RevE.01 30-Nov-2015
                28-Sep-2016
******************************************
%nprocshared=14
Will use up to   14 processors via shared memory.
%mem=43522MB
----------------------------------------------------------------------
#p m062x/genecp opt=(nofreeze,noeigentest,readfc,gdiis,maxstep=10,notr
ust,ts,restart) int=ultrafine freq empiricaldispersion=gd3 5d 7f scrf=
(smd,solvent=acetonitrile)
----------------------------------------------------------------------



......



Berny optimization.
Restoring state from the checkpoint file "/home/gauuser/Gaussian/PhS_Rad_Cyc/4_5
_TS/4_5_TS[Complete_M062X]_Step2[TS_opt].chk".
ONIOM data not found on unit 2.

FileIO operation on non-existent file.
FileIO: IOper= 2 IFilNo(1)=  -997 Len=          20 IPos=           0 Q=  140735655680592


dumping /fiocom/, unit = 1 NFiles =    12 SizExt =    524288 WInBlk =       512
                   defal = T LstWrd =      596480 FType=2 FMxFil=10000

......

dumping /fiocom/, unit = 3 NFiles =     1 SizExt =    524288 WInBlk =       512
                   defal = T LstWrd =       65536 FType=2 FMxFil=10000

Number              0
Base            40448
End             65536
End1            65536
Wr Pntr         40448
Rd Pntr         40448
Length          25088
FileIO operation on non-existent file.
FileIO: IOper= 2 IFilNo(1)=  -997 Len=          20 IPos=           0 Q=  140735655680592
Error termination in NtrErr:
NtrErr Called from FileIO.

按上面方法看似实际有效的报错为“ONIOM data not found on unit 2.”,但其实这根本不是个ONIOM任务,输入文件与ONIOM毫无关系。 其看似报出ONIOM错误的原因是,在正常的opt=restart任务中,程序会在尝试读取 ONIOM data 之后进行标题和Route的读取(见下面的样例),故上述错误很可能是在此处读取这些信息时出错。故也可通过不读取chk文件解决。分析出这个问题需要了解类似的正常任务中“将会有什么输出”,难度较大。

Restoring state from the checkpoint file "/home/gauuser/Gaussian/PhS_Rad_Cyc/2_Z
E_transform/TS[Complete_M062X]_Step2[TS_opt]cont.chk".
ONIOM data not found on unit 2.
Title:  TS [EXTRACT_GEOM]:2,7,12
Route:  #p m062x/genecp opt=(noeigentest,ts,calcfc,nofreeze) int=ultrafine freq
empiricaldispersion=gd3 geom=allcheck 5d 7f scrf=(smd,solvent=acetonitrile) gues
s=tcheck
FncErr=1.00D-07 GrdErr=1.00D-06

还有些杂七杂八的情况,诸如我在做片段初猜的计算时,误将多重度 0 1 0 2 0 -2 输入为 0 1 0 3 0 -3,结果即有如下的错误

(Enter /home/gauuser/g09/l122.exe)
Structure from the checkpoint file:  "sideroidide_CS_1_opt.chk"
FileIO operation on non-existent file.
FileIO: IOper= 2 IFilNo(1)=  -584 Len=          36 IPos=           0 Q=   47012767336112


dumping /fiocom/, unit = 1 NFiles =    33 SizExt =    524288 WInBlk =       512
                   defal = T LstWrd =     1206272 FType=2 FMxFil=10000

Number              0           501           502           503           507           511
Base           216576         43520         72192        180736        181760        196608
End           1206272         44520         76297        180821        181925        197835
End1          1206272         44544         76800        181248        182272        198144
Wr Pntr        216576         43520         72192        180736        181760        196608
Rd Pntr        216576         44520         72192        180736        181760        196608
Length         989696          1000          4105            85           165          1227


......


dumping /fiocom/, unit = 3 NFiles =     1 SizExt =    524288 WInBlk =       512
                   defal = T LstWrd =       65536 FType=2 FMxFil=10000

Number              0
Base            40448
End             65536
End1            65536
Wr Pntr         40448
Rd Pntr         40448
Length          25088
FileIO operation on non-existent file.
FileIO: IOper= 2 IFilNo(1)=  -584 Len=          36 IPos=           0 Q=   47012767336112
Error termination in NtrErr:
NtrErr Called from FileIO.

此时只能根据出错的Link数来推测、测试,得知真正为多重度设定问题。

总结一下我解决这个报错大致思路:

(1)去除所有需要读chk、restart的操作尝试 (2)在大堆数字的上面找到真正的报错,找有效信息 (3)输入文件中用#p代替#,观察出错的Link,思考相关原因 (4)其他杂七杂八问题,硬盘满了,权限问题,路径中有奇怪的字符等等 (5)若在优化过程中出现这个问题,可以将当前结构保存重新开始一个任务,偶尔能解决。 (6)Gaussian Bug,如果排查了上述问题、做了化简,换了新版本、机器、跑多次都能重复这个错误,参阅下文关于Gaussian Bug的说明。

如果做了上述尝试还未能解决这个问题,要向别人提问的话请提供完整的输入输出。

Probably known bugs

▶️ What does a Gaussian Bug look like? ⏳

触发Gaussian的Bug时常有类似这样“模式”的报错:

    Raff turned off since only 57.62% of shell-pairs survive.
    Unable to match L and R vectors in BiOrth.
    Coeff:      0.000D+00 0.000D+00 0.000D+00 0.100D+01
    Logic error in SftOpn.

即以某六字母变量为中心语的某些负面描述 这时就有较大的Bug嫌疑了

另外部分Ntrerr也是Bug

Error termination in NtrErr:
NtrErr Called from FileIO.

遇到“Bug”时

反思自己有没有瞎搞高斯的临时文件。诸如:跑着跑着临时文件被删了;两个文件用同一个CHK、RWF文件名及路径在同时跑(将某个输入文件复制改写为其他结构时常常发生);用脚本比较深层次的控制高斯可能引起的一些奇怪问题(比如用脚本调用 Gaussian external 读梯度时若两个临时文件冲突,会诱发 L601 的 RdWrB1 read garbage pointers 报错。) 首先应更换当前最新Gaussian版本(发稿时为Gaussian09 E.01) 如果仍然报同样问题,可询问高斯客服。 如果不愿询问(如没有版权等),只能更换方法/基组、twitch一下结构、更换积分格点、如果是激发态计算将TD换为TDA等等杂七杂八的方式尝试绕开它了、

Some known bugs:

  • L103,Bad arguments to LAPack or BLAS routine.
  • L401,Diagonalization in DiagDN via DSPEV failed.
  • L906,Internal consistency failure #1 in GetIJB.
  • L914,Unable to match L and R vectors in BiOrth
  • L1002,NIJ > Max2 in MMCore.
Estimated number of processors is:   10
Inverted reduced A of dimension  2103 with in-core refinement.
NIJ > Max2 in MMCore.
Error termination via Lnk1e in l1002.exe at Thu Jun  2 11:12:43 2016.

成因:已知Bug,自G09 C.01版本之后已修复 (G09 C.01 Release Note:“A memory allocation bug for very large systems, which could cause a failure with the message “NIJ > Max2 in MMCore,” was fixed”)

解决:更换G09 C.01或之后版本。

  • L1014,Tx not orthogonal to T
  • Some of the situations involving NtrErr Called from FileIO

Normal outputs that could be mistanken as errors ⏳

Warning!!: The largest alpha MO coefficient is... ⏳

Warning!!: The smallest alpha delta epsilon is... ⏳

(Enter /home/ki/g09/l801.exe)
Windowed orbitals will be sorted by symmetry type.
ExpMin= 2.53D-02 ExpMax= 3.74D+05 ExpMxC= 1.18D+03 IAcc=3 IRadAn= 5 AccDes= 0.00D+00
HarFok: IExCor= 205 AccDes= 0.00D+00 IRadAn= 5 IDoV=-2 UseB2=F ITyADJ=14
ICtDFT= 12500011 ScaDFX= 1.000000 1.000000 1.000000 1.000000
Largest valence mixing into a core orbital is 4.82D-04
Largest core mixing into a valence orbital is 1.44D-04
Range of M.O.s used for correlation: 6 101
NBasis= 101 NAE= 9 NBE= 9 NFC= 5 NFV= 0
NROrb= 96 NOA= 4 NOB= 4 NVA= 92 NVB= 92


**** Warning!!: The largest alpha MO coefficient is 0.21049552D+02
**** Warning!!: The smallest alpha delta epsilon is 0.84363472D-01


Leave Link 801 at Mon Sep 8 14:14:12 2014, MaxMem= 805306368 cpu: 0.7
(Enter /home/ki/g09/l804.exe)

这两组语句是在SCF结束后发现最大分子轨道的系数太大、以及分子轨道间的能量差太小。这一检查的意义是在post-HF(如MP2)、TDDFT等计算中存在以分子轨道系数为被除数、分子轨道间能量差为除数的项,所以分子轨道系数太大、或分子轨道Gap接近0时,这一除法的结果将得到一个较大数值、从而可能引发数值不稳定性,导致结果有误。故而设置这一警告。在基组存在近似的线性相关问题时,如使用了弥散函数,出现这一警告的概率较大。

但这个警告的阈值被设置的较为严格,在正常的计算中也会报出这一警告。如果确认结果基本合理,则可以忽略这两个Warning。 若结果不合理,可:

(1)自己更换或修改基组,尤其是在不需要弥散函数时去掉弥散函数,防止出现接近线性相关的问题 (2)IOp(3/59)可能可以控制检查基组线性相关问题的阈值,可将其放宽尝试 (3) 参见上文中关于“L801,Fatal problem: The smallest alpha delta epsilon is XXXXX”的解决方案。 某些错误常会在该提示后面发生,而Warning前后有空行,颇引人注意,因而常有人误以为这是错误的根源。例如如下段落的报错是“EpsInf not defined for this solvent.”,而不是几个硕大的Warning。

**** Warning!!: The largest alpha MO coefficient is 0.30204407D+03


**** Warning!!: The smallest alpha delta epsilon is 0.88447585D-01


**** Warning!!: The largest beta MO coefficient is 0.30099143D+03


**** Warning!!: The smallest beta delta epsilon is 0.91819897D-01

Leave Link 801 at Fri Nov 6 19:09:46 2015, MaxMem= 1048576000 cpu: 0.3
(Enter /software/gs09/g09/l1002.exe)
Minotr: UHF open shell wavefunction.
NEqPCM: Using equilibrium solvation (IEInf=0, Eps= 13.8200, EpsInf= 0.0000)
EpsInf not defined for this solvent.
Error termination via Lnk1e in /public/software/gauss/g09/l1002.exe at Mon May 11 07:58:05 2015.

附:高斯官方对此问题的回复:

Gaussian tests the values of the MO coefficients and orbital energies after each SCF calculation. When you get a warning about the smallest delta epsilon, the warning is referring to the difference between orbital energies. Since some of the equations in post-HF  methods include a term in the numerator which has the MO coefficients  to some positive integer and a difference in orbital energies in the  denominator, the calculations can become numerically unstable if the  numerator is too large (e.g. MO coefficients too large) or the  denominator is too small (orbital energy difference too small). So, this is a warning that the post-HF results may be affected. However, the criteria for these warnings are fairly strict, so the vast majority of the calculations that have these warning messages do not actually  suffer from numerical instabilities. If you have results that don't seem to make sense, and these warnings are present, it would be a good  idea to repeat the calculation with a different basis set (these  problems tend to be worse for basis sets with many diffuse functions) to verify the results.

Error on total polarization charges ⏳

Cycle  11  Pass 1  IDiag  1:
RMSU=  6.48D-09    CP:  1.00D+00  1.01D+00  1.00D+00  1.05D+00  1.06D+00
                    CP:  9.80D-01  9.86D-01  1.06D+00  1.10D+00  9.43D-01
E= -12097.4433058840     Delta-E=        0.000000000164 Rises=F Damp=F
DIIS: error= 4.95D-08 at cycle  11 NSaved=  11.
NSaved=11 IEnMin=10 EnMin= -12097.4433058841     IErMin=11 ErrMin= 4.95D-08
ErrMax= 4.95D-08  0.00D+00 EMaxC= 1.00D-01 BMatC= 9.84D-12 BMatP= 5.98D-11
IDIUse=1 WtCom= 1.00D+00 WtEn= 0.00D+00
Coeff-Com: -0.204D-04 0.138D-05 0.105D-03 0.650D-03 0.539D-03-0.301D-02
Coeff-Com: -0.154D-01-0.324D-01 0.431D-01 0.322D+00 0.684D+00
Coeff:     -0.204D-04 0.138D-05 0.105D-03 0.650D-03 0.539D-03-0.301D-02
Coeff:     -0.154D-01-0.324D-01 0.431D-01 0.322D+00 0.684D+00
Gap=     0.206 Goal=   None    Shift=    0.000
RMSDP=3.88D-09 MaxDP=5.72D-07 DE= 1.64D-10 OVMax= 1.85D-06

Error on total polarization charges =  0.03783
SCF Done:  E(RM052X) =  -12097.4433059     A.U. after   11 cycles
            NFock= 11  Conv=0.39D-08     -V/T= 2.0073
KE= 1.200978298414D+04 PE=-6.824880392971D+04 EE= 2.436187738767D+04
SMD-CDS (non-electrostatic) energy       (kcal/mol) =      -8.14
(included in total energy above)
Leave Link  502 at Sat Feb 11 02:31:20 2017, MaxMem= 18352963584 cpu:      8769.1

如上面完全正常的任务输出段落所示,这两个语句都不是报错。 认为它是报错其实是英语问题。。。此处 Error 分别指总极化电荷和DIIS的误差,而不是错误。

This type of calculation cannot be archived. ⏳

(Enter /home/gauuser/g09/l9999.exe)

This type of calculation cannot be archived.

Job cpu time:       1 days  8 hours 55 minutes 51.4 seconds.
File lengths (MBytes):  RWF=   1896 Int=      0 D2E=      0 Chk=    138 Scr=      1
Normal termination of Gaussian 09 at Sat Feb 11 03:43:47 2017.

此“报错”出现于任务正常结束前。“Archive” 表示的是Gaussian在一般任务的最后常会输出一大段乱码一样的文本、用于存储任务的各种关键信息。这里“cannot be archived”仅表示 Gaussian 对于 IRC 等类型任务未使其输出存档段信息。

典型的Archieve段落:

1|1|UNPC-LIYUANHE-PC|FOpt|RB3LYP|6-311+G(d,p)|C1H4|LIYUANHE|10-May-201
1|0||# opt freq b3lyp/6-311+g(d,p) geom=connectivity||CH4 OPT FREQ||0,
1|C,6.705202501,1.156069305,0.|H,7.0687325183,0.127811533,-0.000000934
4|H,7.0687519538,1.6701910498,0.8904935408|H,7.0687504279,1.6701921288
,-0.8904935408|H,5.614575104,1.1560825085,0.0000009344||Version=IA32W-
G09RevA.02|State=1-A1|HF=-40.5339325|RMSD=6.153e-009|RMSF=9.920e-005|D
ipole=0.,0.,0.|Quadrupole=0.,0.,0.,0.,0.,0.|PG=TD [O(C1),4C3(H1)]||@

No special actions if energy rises. ⏳

(Enter gpfs/share/home/g16/l502.exe)
Keep R1 ints in memory in canonical form, NReq=136241796.
FoFCou: FMM=F IPFlag=           0 FMFlag=           0 FMFlg1=           0
         NFxFlg=           0 DoJE=F BraDBF=F KetDBF=F FulRan=T
         wScrn=  0.000000 ICntrl=       600 IOpCl=  0 I1Cent=           0 NGrid=           0
         NMat0=    1 NMatS0=  14878 NMatT0=    0 NMatD0=    1 NMtDS0=    0 NMtDT0=    0
Symmetry not used in FoFCou.
Two-electron integral symmetry not used.
Closed shell SCF:
Using DIIS extrapolation, IDIIS=  1040.
NGot=  9895542784 LenX=  9784774425 LenY=  9784735568
Requested convergence on RMS density matrix=1.00D-08 within 128 cycles.
Requested convergence on MAX density matrix=1.00D-06.
Requested convergence on             energy=1.00D-06.
【No special actions if energy rises.】

Cycle   1  Pass 1  IDiag  1:
Inv3:  Mode=1 IEnd=     1030188.
Iteration    1 A*A^-1 deviation from unit magnitude is 3.00D-15 for    375.
Iteration    1 A*A^-1 deviation from orthogonality  is 2.66D-15 for    435    266.
Iteration    1 A^-1*A deviation from unit magnitude is 2.66D-15 for    375.
Iteration    1 A^-1*A deviation from orthogonality  is 2.66D-15 for    379    156.

此语句不是报错,它只是在描述SCF过程中遇到该步骤的能量比上一步骤高的时候,程序该做什么。 这一语句的选项可以由 IOp(5/86) 控制,例如设置 IOp(5/86=101202) 时会输出下列内容(默认为 IOp(5/86=101100)):

Reduce DIIS space if energy rises from previous iteration.
Dynamic level shift is off after energy rises.

初学者常以为这句话是报错的原因是程序在输出这句话之后会进行SCF迭代,耗时较长,如果初学者不加#p输出SCF迭代的详细信息,那输出文件末尾会保留在“No special actions if energy rises”很久,令初学者误以为出错了。平时计算时都应该加上#p来输出具体SCF信息。

DIIS: error ⏳

End of XXXXXX F.D. properties file XXX does not exist. ⏳

FullF1:  Do perturbations      1 to       3.
Isotropic polarizability for W=    0.000000      466.35 Bohr**3.
End of Minotr F.D. properties file   721 does not exist.
End of Minotr F.D. properties file   722 does not exist.
End of Minotr F.D. properties file   788 does not exist.
Leave Link 1002 at Fri May  3 18:59:51 2019, MaxMem= 10995105792 cpu:          153536.8 elap:            4809.4
(Enter /gpfs/share/home/1501110295/g16/l601.exe)
G2PCM: DoFxE=T DoFxN=T DoGrad=T DoDP/DQ/DG/TGxP=FFFF NFrqRd=   0 IEInf=0 SqF1=F DoCFld=F IF1Alg=4.
GePol: Maximum number of non-zero 1st derivatives   =     301
End of G2Drv F.D. properties file   721 does not exist.
End of G2Drv F.D. properties file   722 does not exist.
End of G2Drv F.D. properties file   788 does not exist.
Leave Link 1110 at Fri May  3 06:07:58 2019, MaxMem= 10995105792 cpu:           17529.5 elap:             549.1
(Enter /gpfs/share/home/1501110295/g16/l1002.exe)

这几个语句的含义不明,但不是报错。其常被用户误以为是报错的原因是其出现在计算频率时较为耗时的进程前面,对较耗时的体系没有新的输出,故而使用户误以为“跑到这里就卡死了”,但其实经常是高斯正在进行计算,确认高斯还在运行中之后等待算完即可。

Abnormal situations that doesn't give an error message

Abnormal situations that don't give an error message

▶️ IRC job stopped prematurely, or goes in the other direction ⏳

image

▶️ One or more small imaginary frequencies after optimization ⏳

有时优化后会出现一小虚频,虚频数值在几至几十cm^-1。 image 解决:

若满足:

(1) 使用的是Gaussian 09 E.01之前的版本,如Gaussian 09 D.01 (2) 使用了D3BJ色散校正(关键词为empiricaldispersion=GD3BJ)

则可选择在优化和频率计算中同时: (1) 换用Gaussian 09 E.01或之后版本 (2) 换用empiricaldispersion=GD3 (3) 不用色散校正

这个Bug的具体解释见:http://bbs.keinsci.com/thread-2410-1-1.html

若不满足上述两个情况:

看《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》(http://sobereva.com/278)

▶️ G09 Windows exit at L1.exe without an error message ⏳

image image 成因:目前网上能找到的G09W都是32位程序,故最多可设置使用约1300MB内存(%mem=1300MB,实际极限在1380~1600之间),4 核并行,16GB 硬盘(切割Scratch file无效)。如果超过此值可能触发不报错退出。如上图即为故意将内存设为 %mem=1500MB 时的“无报错错误输出”。

解决:

将内存、核数、硬盘设置到相应极限以下 有钱的去买64位的Windoes G09 程序( ╮(╯▽╰)╭ 然后打个包给大家就更好了) 换用Linux下的Gaussian程序

▶️ Abnormal bonding / unbonding in GaussView ⏳

成因:首先“键”并不是原生的量化概念,量化的输入、输出也(几乎)不依赖键连信息,显示成键情况只是为了辅助理解而已。很多非经典的情况下经典意义上的“键”并不是一个好的概念。

      Gaussview 中的显示的成键状态(在未人为规定时)是用原子间距离判断的,即两原子间键长在这对儿元素(比如 C 和 C )对应的某个范围内(如145 pm-165 pm),即判定为相应键级的成键(如C-C单键),这个成键情况的判断显然是十分粗糙的,不必在意。

常见的成键、断键异常有这几种情况 下图中H2O2是优化后,O-O 键键长变得较长,显示 O-O 之间断键了,但测量 O-O 原子间距离可知仅 145pm,做进一步的其他键级分析也可以得出其中显然有一个共价键,故 “O-O断键” 只是显示问题而已,并非结构优化等过程有错。 中间的例子是一个芳环,芳环键的判定常有错误(不显示为共轭的“1.5级键”),经常判断为全是双键、或完全混乱等等,同样不必在意。 第三个例子较为特殊,其打开的是一个fchk(或chk)文件。在优化过程中,Gaussian会在fchk文件中保存第一步优化/IRC 时判断出的成键信息 / 或 geom=connectivity 定义的连接信息,故不论后续优化过程如何进行,打开时GaussView都会显示优化/IRC 刚刚开始时的成键信息。第三个例子就是一个IRC过程,在初始结构中,1-2两个氮原子很近,根据原子间距离判断图中1-2两个氮原子间有单键,IRC/优化后虽然 1-2 两个氮原子远离,但第一步的键连关系被保留,显示为一根异常长的单键。观察这个结果中的键长、键角,根据常识判断是合理的,说明其也是正确的,不说明计算有问题。如果打开的是Log文件,且没有记录geom=connectivity则没有这个显示问题 image 解决:GaussView (及其他很多软件,如Spartan)显示的成键情况很粗糙,不论是否出现未预期的结果,都应自己通过化学常识判断几何优化过程、显示的成键情况等是否合理。另外可以用Multiwfn计算Mayer键级来辅助判断。

对于上图第三种情况,可以使用下图所示的Rebond工具,让GaussView根据当前键长重新粗略判断成键情况,可以去除一些明显错误的异常成键 image

▶️ Viewing IRC jobs in GaussView, it gives an abnormally low energy point at the X=0 point ⏳

症状为在 GaussView 中打开 Result - IRC/Path 时,在反应坐标零点处附近出现一低能点,使图形看起来是“第一个点就掉下去了”: image

成因:GaussView 显示 Bug。

      仔细观察可发现,如下图所示,只要将“掉下去的”点移动到曲线的一端(红叉所示,左右不定),曲线就正常了(绿线所示)。
      ![image](https://user-images.githubusercontent.com/18537705/160455534-05b91e6d-0355-45ce-891b-4d403ec699ac.png)

其原因在于 GaussView 在读取未正常结束的(尚未跑完、跑出错)的 IRC 时,当前点的反应坐标尚未计算出(如下面段落所示的信息缺失),但当前结构及其能量却已经得到并被GaussView读取,故GaussView 默认其横坐标为0,导致点的顺序出错,显示为这一问题。

     CHANGE IN THE REACTION COORDINATE =    0.30843
   NET REACTION COORDINATE UP TO THIS POINT =    2.77474
  # OF POINTS ALONG THE PATH =  19
  # OF STEPS =   1

解决:如果任务正在运行,使其跑完 如果任务已经跑完,应是出了什么问题没有正常结束,看输出文件具体报错解决。

▶️ "Missing or bad data: Alpha Orbital Energies" when opening chk or fchk files with GaussView ⏳

解决:两种方法 (1) 先将chk转换为fch文件,用Multiwfn软件载入fch文件,用主功能100-选项2转换为fch文件给GaussView读入。或者:(2)先将chk转换为fch文件,然后用文本编辑器打开fch文件,看看里面的Number of basis functions和Number of independent functions是否一致。如果不同,将independent改成independant

Other common entry-level questions

▶️ “How long will this calculation takes?” ⏳

常有人给出一个体系询问这个计算大概需要多长时间,这种问题往往难以回答。

计算所需的时间取决于很多因素,包括而不限于下列因素(请注意,只是举例,即使依次回答下面的问题也常不足以给出靠谱的猜测): 体系的大小(原子数与元素构成) 方法(除post-HF和DFT/HF这样的比较外,还例如开RI时杂化泛函比纯泛函慢) 基组数量、种类(相同数量的弥散函数远比极化函数耗时多);轨道数 电子态及电子结构的复杂程度(例如一般分子需要10~20步SCF,稀土有时需要几百步SCF;同样的基组数,一般有机物也远比金属簇容易) 任务类型 初猜的好坏 任务的顺利程度(如优化步数、SCF步数) 加速或“减速”设定(如RI, int=ultrafine等); 电脑的CPU(核数、频率、架构) Gaussian版本(例如avx2就比avx版本快一些,更极端的还可能涉及数值导数和解析导数的差别) 内存、硬盘情况(大小、速度)及为任务分配了多少,是否有其他任务在共享,使用N个核的超线程时是用了N个线程还是2N个线程等等等等琐碎的问题。

以上任何一个因素改变都可能显著的、在数量级层面上改变计算所需的时间,同时也不会有人了解所有的机器环境(用过各种CPU之类)、所有的任务和体系。故往往此类问题根本没法回答。

比较合理的估计耗时的方法是在自己机器上测试计算中某个步骤的耗时,用它估算。有限的可供参考的“感觉”大致如下: (我主要用(非双杂化的)DFT,如B3LYP、M06-2X之类,下面仅对这类方法、仅对一般有机体系适用。阅读下文时请注意区分“一步SCF迭代”,“(整个)SCF过程”,“一步几何优化”,“整个几何优化过程”的区别)

一次单点的耗时,可以通过一步SCF迭代进行估计。初次SCF一般在十几步至几十步之间收敛,超过50步应开始检查SCF震荡;优化时后续优化步骤的SCF会读取之前的初猜,故优化的第二步开始,SCF所需的步数常会显著减少,一般在几次至十几次(常可见优化中第一步L502耗时最多,后续L502耗时较少)。可根据这一经验、通过一步SCF迭代大致乘出整个SCF过程的耗时。 解析梯度计算的耗时很少,故几何优化的耗时(等于一次单点+一次梯度,不含calcfc)一般只比单点高一点,可略。 几何优化任务,对小体系一般在二三十步以内优化收敛,大体系(100+原子)在50~150步都有可能,再多需谨慎考虑优化震荡的可能性。另外上面说过优化中第一次SCF过程常常比后续SCF所需的时间长很多,所以估计几何优化的耗时使用第二步几何优化的时间,乘以可能需要的几何优化步数来估计比较合适。 频率计算一般需要单点几倍至几十倍的时间,少则3倍,多则近100倍也遇到过,这一倍数似乎与几何的变量数目有关系。可用这个规律估计各处计算频率的耗时(含优化完成后、使用calcfc时优化第一步、以及calcall的每一步涉及的频率计算)

其他可能(严重)影响耗时的情况:

SCF、优化震荡:做大任务时,应设法实时监测SCF和优化的步骤,可以自己写一些小工具,不要提交个任务,跑了两个星期跑到500步报错了,才发现从20步开始就在震荡了。 其他理论方法:一些有解析梯度、频率的 post-HF方法(含双杂化泛函如B2PLYP),其解析梯度和频率耗时与单点的比例一般高于HF和DFT方法,解析梯度可至单点的几倍,解析频率 复杂电子结构:有的体系的电子结构复杂,诸如含镧系原子的分子,其SCF收敛可能需要几百步,远多于一般体系,应特殊处理(此时优化的时候用calcall甚至可能比不用的时间还短) 使用特殊选项:诸如使用 scf=xqc时,有时可以用常规方法收敛,有时不行;如果使用了qc,会使得一次SCF收敛的时间大大增长,以往的“SCF需要多少步”的经验也不再适用 数值梯度和数值频率:对于没有数值梯度的方法,诸如CCSD(T),每步几何优化的耗时等于6N步单点能(N为原子数),此时梯度计算所需的时间显然不可忽略;对没有解析Hessian的方法,计算一次频率的耗时为6N次[单点+梯度]

对这个问题有不同或其他意见欢迎补充

▶️ How to generate an initial guess for transition state geometry optimization with flexible scan ⏳

该方法是寻找过渡态的备用方法,并非一定成功、且耗时甚高,这绝不是找过渡态的“标准方法”,使用者稍有提供初猜的经验时、绝大部分过渡态都可以用opt=TS寻找到。

柔性扫描是有效的获得过渡态初猜的傻瓜方法,在初学者刚接触某类反应、不了解大致的过渡态结构(而且还懒得去查文献的时候)时、用此法不必动很多脑筋、比较robust;在研究有些较难找到的过渡态,比如某些自由基过程时也比较有效。柔性扫描很简单、但却常有人问,故在此说明。请务必注意本段最后的注意事项。

我们以下列自由基氢转移反应为例:

image

在GaussView中建模,首先创建一个甲烷,使用Add-Valence按钮点击其中一个H,将新出现的氢原子取代为SH基团。基本建模就完成了。

image

随后需要调整键长。这个反应中涉及两个键、分别是 S-H 和 H-C 键,但我们只需扫描其一,因为在优化过程中另一个键长会随之弛豫变化。

我个人建议,对新手来说、如果自己对过渡态中的键长没有较好把握的话,做两次扫描,均以自己猜想的过渡态为起点,这有两点好处,一是可以在扫描过程中监测扫描过程,发现已有最高点时可立即停止计算,节约计算量;二是这样在优化过程中另一边不会“飞掉”,如本例中,若 C-H 键键长较小,这个体系就变成了分子间弱相互作用体系,只用一次扫描时容易在开始阶段优化不收敛、或因为初猜原因跑到不希望要的构型上。

这个反应中我们以扫描反应中涉及的 C-H 键为例,正常得C-H键键长在110 pm左右,过渡态中势必会较长,假设我们猜过渡态中 C-H 键键长为 1.7A 左右。

image

随后在 Edit-Redundant Coordinate 中设置柔性扫描 C-H 键,扫 5 步,每步 15 pm (即 170 ~ 245 pm 的 6 个结构),用 B3LYP/6-31G(d) 为例提交运算;再重新设置每步为 -15pm,扫描 3 步(即 125 ~ 170 pm 的 4 个结构),提交运算。

image

GaussView产生的第二个输出文件如下:

%nprocshared=2
%mem=1200MB
%chk=C:\Users\LiYuanhe\Desktop\1.chk
#p opt=modredundant freq b3lyp/6-31g(d)

Title Card Required

0 2
C                 -3.01654961   -1.31556721   -0.25719647
H                 -2.65989518   -2.32437721   -0.25719647
H                 -2.65987677   -0.81116902   -1.13084798
H                 -4.08654961   -1.31555402   -0.25719647
H                 -2.44987313   -0.51418690    1.13084798
S                 -2.01319891    0.10334733    2.20045870
H                 -2.44953630    1.33854341    2.20102655

B 5 1 S 5 0.150000

其中 opt=modredundant 是进行柔性扫描的关键词,分子段落后面的 “B 5 1 S 5 -0.150000” 表示增加一个冗余内坐标,内容是一根化学键(B),在 原子 5 和 原子 1 之间,模式为扫描(S),共扫描 5 步,步长为向减小键长方向的 0.15 A。

计算完后可在 GaussView 中 Result-Scan 分别查看两次任务的曲线。为方便考察,这里将两个结果做到一条曲线上:

image

可见其最高点在 1.7 A 处。故我们在 GaussView 中取 1.7 A 的扫描结果作为过渡态初猜,用opt=TS,在同一水平下找过渡态,近4步优化就找到了正确的过渡态,过渡态 C-H 键键长 1.675 A。

%nprocshared=2
%mem=1200MB
#p opt=(calcfc,ts,noeigentest) freq b3lyp/6-31g(d)

Title Card Required

0 2
柔性扫描最高点的分子坐标

(如果你完全按上述设置,会发现扫描时最左侧一点时会遇到了L9999 的 opt 未收敛错误,可以按照本帖上面提到的方法解决,但也可无视,因为此时打开输出文件可见路径上的最高点已经找到了,就不用算更多点了)

几点说明: 该方法应作为寻找过渡态的备用方法(比如上面反应的自己画初猜用 opt=TS 其实也容易找到),而不是每次找过渡态都用它。该方法较为耗时,但还可以接受。例如扫描时取10个点时,因为后续的优化是用前面的结构做初猜的,其并不等于做十个优化任务,而只有第一步是完整的优化、后续优化步大概只需要一般优化任务的一半甚至更少的步数和时间,故柔性扫描10个点一般也只是25倍优化任务所需的时间,对不十分大的任务还是容易接受的。 需强调使用的是柔性扫描(opt=ModRedundant),不是刚性扫描(Scan),后者对找过渡态几乎没用。 做柔性扫描的级别可以比实际寻找过渡态的级别稍低,但不能太差。例如过渡态用PBE0/def-TZVP计算,则柔性扫描最低可以在 PBE0/def2-SV(P)、或B3LYP/6-31G(d)级别下做;而如果用HF/6-31G(d)这样的级别耗时不低、得到的结果却帮助有限。 柔性扫描只是帮助产生过渡态初猜,而不是“寻找过渡态”。柔性扫描的最高点(即使步长无限小)也几乎都不是过渡态。柔性扫描的曲线有时与IRC相差甚远。 步长 10 pm 左右基本是够用的,如果势能面比较复杂,对初猜十分敏感,也可以先用较宽的步长扫描,然后对感兴趣的区间用更小的步长扫描(如上述例子中、可以最后再用0.03的步长在1.51.8之间扫描) 本文的扫描和过渡态寻找都用 B3LYP/6-31G(d) 只是举个例子,实际计算中应当选取合适的计算级别。扫描和寻找过渡态可以用不同的级别,但如果两种方法的势能面差别太大,可能效果不好,若体系、级别很耗时,可如先用 B3LYP/6-31G(d) 做上述扫描,最后再用高级别在最高点附近确认三个点,如在本例中用 M06-2X/def-TZVP 重新确认 1.6 A, 1.7 A, 1.8 A 三个点,确认大小关系正确(方法是在 Edit-Redundant Coordinate 中选择 Freeze 而不是 Scan)。 注意柔性扫描过程中会出现因结构跳变造成的Artifact,体现为曲线不平滑、有“悬崖”,如下图所示。注意此时曲线上的极大值点与过渡态毫无关系,需为连续、平滑的曲线(如上面的例子中所示),才体现过程中确实存在过渡态、且可作为初猜。遇到这种扫描曲线不平滑的情况应视扫描轨迹的具体情况解决问题,尝试避免跳变的发生。

image

▶️ “How to put a charge / electron on this atom / fragment?" ⏳

http://sobereva.com/82 一文第六节

▶️ How to calculate activation-strain (distortion-interaction) energy ⏳

Activation Strain 模型和 Distortion-Interaction 完全是一个东西。前者由 Bickelhaupt (1999) 基于前人的一些想法建立起来,提出 Activation Strain 模型的名称;后来 K. N. Houk 在 2007 年 JACS 上发表了一篇文章,未引用 Bickelhaupt,只引用之前更早的一些文章,称之为 Distortion-Interaction,内容与 Bickelhaupt 是完全相同的,没有任何新东西。但或许只是因为他有名,故而大家多称之为 Distortion-Interaction Model。在此鄙视这种行为。

Distortion-Interaction 模型是一种对活化能进行能量分解的方案。对于一个双/多分子反应(单分子反应不适用),distortion能量定义为过渡态各组分在过渡态构型下、无相互作用时的能量与稳定底物之间的电子能差,此值一般为正值;interaction 定义为前述无相互作用体系与真实的过渡态电子能差,此值一般为负值。二者之和应精确等于活化能。 以下述Diels-Alder反应为例。

image

定义如下能量。E1、E2为相应物质优化后(opt)的电子能。E5为过渡态优化后的电子能(opt=TS)。随后将过渡态优化后的结构提取出来,删去其中的蓝色部分的几个原子,不优化,在同样级别下计算单点能(sp),得到能量E3;类似的删除红色部分,计算单点能得到E4。

image

则有: 活化能 Ea = E5-E1-E2;

Distortion Energy = E3 + E4 - E1 - E2;

Interaction Energy = E5 - E3 - E4。

可见 Ea = Interaction Energy + Distortion Energy,是一个精确地能量分解。

这个能量分解方案简单,不需要任何特殊技术,耗时相对优化也很低。有时可以为分析提供一定角度。

除了对过渡态做这件事之外,也可以对IRC上的每个点做完全一样的能量分解。观察 distortion 和 interaction 在“反应过程中”的变化。

About

Reliable explanation and ways to solve the error messages for the Gaussian quantum chemistry program. This is a per request English translation for a long blog in Chinese I wrote a few years back.

Resources

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published