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