follow up to MOLPRO ROHF convergence difficulties (optimizing a geometry using procedures)

The end result of the last post was that I found level shifts of -0.3, 0.3 worked for my system. A brief tangent for those of you who have perused the manual and saw their suggestion of -1.0,-0.5, I did try that one as well and it does not work in this particular case. Ok, back to the question at hand, how do I optimize a geometry when ROHF won't even converge for a structure that is a local minimum? One way around this is to use MCSCF with one open orbital and run the converged orbitals through ROHF using only one iteration. However, that is a story for another day, in this case I will use the level shifts I listed above. My basic thought was that I really only want to use the level shifts if they are needed, otherwise I will use the default values. These two ideas lead to using a procedure and boolean logic within said procedure.

basis=aug-cc-pvdz

proc autoshift
   rhf
   if(status.lt.0) then
     status,rhf,clear                          !make sure calc does not abort when stock rohf does not converge
     {rhf,shifta=-0.3,shiftb=0.3;        !shifts found using grid search
     start,atden}                              !start with atomic densities again, otherwise rohf will use same "bad" orbitals
  endif
  uccsd(t)
endproc

{optg,proc=autoshift,root=1,energy=1.d-11,step=1.d-5,gradient=1.d-5,displace=symm,numhess=2,hessproc=b3lyphess}

proc b3lyphess
  rhf
  if(status.lt.0) then
    status,rhf,clear                           !same as for autoshift
    {rhf,shifta=-0.3,shiftb=0.3;                                                   
    start,atden}
  endif
  ks,b3lyp
endproc

Here I have made two procedures, one for the optimization using CCSD(T) and one for the Hessian (MOLPRO does far better with frequent numerical hessians when optimizing open-shell doublet molecules) using B3LYP. I will talk about the first one because they are almost identical. Beginning with the if statement, I hate it check if the stock ROHF calculation finished properly. If it did not, then I run ROHF with the level shifts I found in the previous post, taking care to start with the atomic densities and not the orbitals from the previous ROHF calculation. I could put in several of these if statements with various values for the level shifts to be on the safe side, but I wanted to keep this example on the simple side. If this ROHF calculation fails then I still want MOLPRO to error out, so I do not have a status,rhf,clear statement in the if structure. After that I run the UCCSD(T) calculation. The only other bit to include is telling optg that I want to use a procedure, which is accomplished with the proc=autoshift. This is currently running, but it made it past the initial calculation which is progress.

automation of MOLPRO and ROHF convergence testing

Every now and again I run across cases where open-shell systems cause ROHF in MOLPRO to choke. Usually this is easily solved by starting out with STO-3G, stepping up to 6-31G(d), cc-pVDZ, and, finally, aug-cc-pVXZ. MOLPRO seamlessly reads the orbitals from the previous calculation, if available, and projects these onto the new basis set. However, this approach does not work on some occasions and neither does increasing the maximum number of iterations (I have never run into an example where increasing maxit past 100 has helped). This left me with wanting to try level shifts, which would be tedious for anything more than 3 or 4 tries. This leads me to use BASH scripting to automatically generate input files covering a range of level shifts.

for SHIFTA in -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 ; do for SHIFTB in -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 ; do echo $SHIFTA $SHIFTB ; cp testt-a.inp test-shift${SHIFTA}${SHIFTB}.inp ; sed -i "s/shifta=-0.5,shiftb=-0.5/shifta=${SHIFTA},shiftb=${SHIFTB}/g" test-shift${SHIFTA}${SHIFTB}.inp ; done ; done

The nested for loops generate 81 input files, each having a different set of alpha and beta shifts, ranging from -0.5 to +0.3. The echo line lets you know how far along in the process it is. The cp line creates the input file from a template, test-a.inp. Now to the most interesting line where sed replaces the shift values from the template with the values held in variables shifta and shiftb.

The next step is to run all of these calculations. I opted for simplicity and used the same nested for loops.

module load molpro/2012.1.8 ; for SHIFTA in -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 ; do for SHIFTB in -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 ; do echo $SHIFTA $SHIFTB ; molpro -n1 <./test-shift${SHIFTA}${SHIFTB}.inp >./test-shift${SHIFTA}${SHIFTB}.out; done ; done

Finally, one needs to look at the results. I echo the shift values and then the ROHF energy. The gawk bit just gets rid of the grep term I use to pull this energy from the output file.

for SHIFTA in -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 ; do for SHIFTB in -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 ; do echo $SHIFTA $SHIFTB $( grep "RHF STATE 1.1 Energy" ./test-shift${SHIFTA}${SHIFTB}.out | gawk ' { print $5 } ') ; $( grep "No convergence" ./test-shift${shifta}${shiftb}.out ) ; done ; done >./631gs-shift.txt

From here it is easier to change other aspects with one line, such as the basis set, maxit, etc. If I am trying to get the ROHF calculation to converge I should really have grep find convergence issues. On top of that I would like to know which parameters, if any, worked and which were bogus. I modified the above line to add in a grep for "No convergence" and I dump all of the results to a file, which allows me to use the BASH sort function to parse the data either by energy or by convergence. The below line sorts based upon the energy, but one could change 3,3 to 4,4 and sort by convergence.

sort -t$' ' -k 3,3 631gs-shift.txt