State-specific gradients for state-averaged MCSCF and MCSCF state tracking.

 

State-specific gradients for state-averaged MCSCF

 

The approach used by Firefly to compute state-specific (SS) gradients for state-averaged (SA) MCSCF is described in details here and is based on the differentiation of effective gradient vector computed with state-averaged density matrices over the weight of the state of interest. The differentiation is performed using second-order finite differences. Therefore, due to their seminumeric nature, the computed gradients are more sensitive to numerical errors in the computed solution of SA-MCSCF than the ordinary MCSCF gradients are. Hence, one needs to increase the overall precision of computations as will be discussed below. Computationally, state-specific gradients for SA-MCSCF are approximately two to thee times more costly than gradients for MCSCF without state averaging. The MCSCF procedure itself can be arbitrary i.e. it can be of any supported type, can use both GUGA and ALDET CI code, and can utilize any of available MCSCF convergers. However, the use of Firefly’s unique SOSCF converger is strongly recommended. The relevant input groups are $mcscf and $mcaver

 

The $mcscf input group contains several relevant keywords, namely istate, ntrack, acurcy, and engtol.

istate                    the state number for which SS gradient will be computed. States are numbered starting from one. For example, $mcscf istate=2 $end will compute gradients for the first excited state of given multiplicity and symmetry type as specified by other sections of input file. The only exception is the case of ALDET CI with pures option disabled, for which the numbering will include all states regardless of their multiplicity. One should not use a nonzero istate option for MCSCF runs without state-averaging or when state-specific gradients are not required! The default value is istate=0 i.e. state-specific gradients are disabled.

ntrack                  if nonzero, activates MCSCF states tracking and remapping feature of Firefly. More precisely, ntrack defines the number of lowest roots (states) to be tracked and, if necessary, remapped to other states. The present implementation of state tracking in Firefly is rather simple and is based on the analysis of overlap matrix of current states with the previously computed ones.  See the description of $track group below for more details and for additional keywords.  The default for ntrack is ntrack=0 i.e. do not track states at all. State tracking is cheap and can be of great help e.g. during geometry optimization of excited states as it can detect root flipping and remap states accordingly so it is generally a good idea to activate it.

acurcy                   the major convergence criterion for MCSCF,  the maximum permissible asymmetry in the Lagrangian matrix. While its default value is 1.0E-05 and is suitable for single-state MCSCF, computation of SS gradients for SA-MCSCF requires tighter convergence.  The recommended values of acurcy are in the range 1.0E-07 - 1.0E-08.

engtol                  the  secondary convergence criterion, the MCSCF is  considered converged when the energy change is  smaller than this value. The default is engtol=1.0E-10. Again, SS gradients require tighter convergence; the recommended values of engtol are in the range 1.0E-12 - 1.0E-13.

 

The $mcaver input group defines the details on how exactly SS gradients are computed.  The relevant keywords are deltaw, conic, and hpgrad. Besides these, this group contains additional keywords controlling parts of Firefly’s code used for location of interstate crossings and conical intersections, namely:  ssgrad, jstate, a, b, target, xhess, xgrad, multiw, shift, and penlty. These keywords will be described elsewhere.

deltaw                 the step size (dimensionless) over state’s weight used in finite differencing. The recommended values are in the range from 0.0005 to 0.0025. The default value is deltaw=0.0015.

conic                     selects one of the three programmed approaches to finite differencing.  Valid values are conic=0, 1, or 2. The default is conic=0

Conic=0 selects the use of central (i.e. symmetric) second order finite differences.  The calculations are performed as follows. First, calculations on SA-MCSCF energies and effective gradient are performed with the original weight of state # istate increased by deltaw.  Second, calculations on SA-MCSCF energies and effective gradient are performed with the original weight of state # istate decreased by deltaw. Third, the calculations on SA-MCSCF energies are performed using unmodified original weights, and finally, state-specific expectation value type density matrix is computed for state # istate.  This is the most economical way of computations.  In addition, it provides the way to obtain the state-specific properties for the state of interest.

Conic=1 selects the alternative approach based on the use of forward finite differencing scheme. The latter is more stable in the case of nearly quasi-degenerated CI roots and hence it is more suitable for location of conical intersections. The calculations are organized as follows.  First, calculations on SA-MCSCF energies and effective gradient are performed with the original weight of state # istate increased by deltaw. Second, calculations on SA-MCSCF energies and effective gradient are performed with the original weight of state # istate increased by deltaw/2.  Finally, the calculations on SA-MCSCF energies and effective gradient are performed using unmodified original weights but state-specific density matrix is not computed hence no state-specific properties are available.   

Finally, conic=2 selects another approach which is similar to conic=1 but is even more robust in the vicinities of conical intersections.  First, calculations on SA-MCSCF energies and effective gradient are performed with the original weight of state # istate. Second, calculations on SA-MCSCF energies and effective gradient are performed with the original weight of state # istate increased by deltaw/2.  Finally, the calculations on SA-MCSCF energies and effective gradient are performed with the original weight of target state increased by deltaw. Similarly to conic=1 no state-specific density matrix is computed thus the state-specific properties are not available.

hpgrad                logical variable. If set, requests extra high precision during computations of the two-electron contributions to the effective gradients.  This may significantly slow down computations and usually does not considerably increase the precision of the computed state-specific gradients.  This is why this option is disabled by default (hpgrad=.false.).

 

Hints on increasing the overall precision of SA-MCSCF calculations

 

Typically, the separate stages involved into MSCSF procedure are: evaluation of two-electron integrals, integral transformation, CI procedure, computation of one- and two-particle density matrices, and the orbital improvement step. For SS gradients, one needs to increase the precision of each of these steps. Here we provide a synopsis of related input groups and keywords and give some recommendations on the optimal values for the latter.

1.       Precision of two-electron integrals is controlled by inttyp, icut, and itol keywords of the $contrl group as described in the manual.  For SS gradients, the recommended values are inttyp=hondo, icut= at least 11 (value of 12 or 13 is even better), and itol=20

2.       Precision of integral transformation stage is controlled by the cuttrf keyword of $trans group as described in manual.  For SS gradients, the recommended value of cuttrf is 1.0d-13 or tighter.

3.       Precision of CI step is controlled by a single keyword for ALDET CI code and by two keywords for GUGA CI code. For ALDET code, it is cvgtol of $det group. The recommended value of cvgtol is 1.0d-7 to 1.0d-8 to 1.0d-9. For GUGA CI, these are cutoff of $gugem group and cvgtol of $gugdia group. The recommended values are cutoff=1.0d-20 and cvgtol=1.0d-7 to 1.0d-8 to 1.0d-9.

4.       For GUGA CI, precision of two-particle density matrix computation step is controlled by cutoff keyword of $gugdm2 input group. The recommended value is cutoff=1.0d-15 or tighter. For ALDET CI, there are no additional keywords required.  

5.       The precision of the orbital improvement step can be improved using several keywords of $moorth and $system groups.  It is recommended to set:   

$system kdiag=0 nojac=1 $end

$moorth

nostf=.t. nozero=.t. syms=.t. symden=.t. symvec=.t. symvx=.t.

 tole=0.0d0 tolz=0.0d0

$end

 

An example

 

An archive with commented sample input and output files for both single-state MCSCF gradient and SS gradient for SA-MCSCF computations can be found here

 

MCSCF state tracking.

 

The present implementation of MCSCF state tracking in Firefly is based on the analysis of overlap matrix of the current MCSCF states with the previously computed and possibly already remapped states which serve as the reference vectors.  First, the overlap matrix is computed every MSCSF iteration. The diagonal elements of this matrix are then scaled to decrease the probability of false root flipping detection in the regions where MCSCF states undergo rapid changes i.e. near avoided crossings or conical intersections. Finally, the modified overlap matrix is decomposed using a dedicated procedure and the optimal new mapping for the new states is computed. This scheme is very cheap and simple and typically works well but cannot ideally handle all possible situations. It may be modified or retuned in the future Firefly versions.  Note, tracking will abort the job if the number of available CI vectors is less than the number of MCSCF states to track. The latter is defined by the $mcscf ntrack= option as described above.

The $track input group is used to control the details of MCSCF state tracking. This group contains the following keywords:

tol                         The scaling factor for diagonal of the overlap matrix. The default value is 1.2

update                Logical variable. If set, reference vectors will be updated and replaced by the remapped current CI vectors at the end of each MCSCF iteration. If not set, the CI vectors from the very first MCSCF iteration at the initial geometry will be used as the reference vectors throughout all calculations. The default is not to update reference vectors i.e. update=.false.

freeze                Logical variable. If set, the final remapping scheme of the first of three MCSCF calculations used to compute SS gradient for SA-MCSCF will be applied “as is”  during the second and third stage of gradient computations.  Otherwise, the dynamic tracking will be active throughout all three MCSCF procedures.  Normally, this flag must be set to get reliable results is the MCSCF states are quasi-degenerated. The default is freeze=.true.

reset                  Logical variable. If set, state tracking will be reset at the beginning of each MCSCF computations and the reference vectors will be re-initialized by the vectors from the first CI step. Default is reset=.false.

sticky                 Logical variable to be used together with reset option. If reset and sticky are both set, state tracking will be reset at the beginning of each MCSCF computations but the existing reference vectors will be reused. Default is sticky=.false.

delciv                  Logical variable. At the end of each CI procedure, converged CI vectors are stored in the special file. If delciv is set to .true., the file with converged CI vectors of previous CI stage will be deleted before each new CI step so that these old vectors  will not be used as the initial guess for the new CI procedure. The default value is to reuse old CI vectors as the initial guess and thus to keep the file intact i.e. delciv=.false.