Constrained geometry optimization does not converge

Questions
#1

I want to perform a constrained geometry optimization job with an ethanol+water system, with the ethanol molecule fixed and the water geometry (bond lengths and angle fixed), so as to only optimize the ethanol water distance and rotation of water molecule. This is the input I am using,

$rem
JOBTYPE=OPT
METHOD=HF
BASIS=6-31G*
SCF_CONVERGENCE=8
GEOM_OPT_DRIVER=OPTIMIZE ! Use the old driver
GEOM_OPT_MAX_CYCLES=200
NO_REORIENT=TRUE
SYM_IGNORE=TRUE
HARM_OPT=1 ! Turn on harmonic confining potential
HOATOMS=9 ! Number of confined atoms
HARM_FORCE=450 ! Force constant of the potential
$end

$harmonic_opt
1 2 3 4 5 6 7 8 9
$end

$opt
CONSTRAINT
stre 10 11 0.9572
stre 10 12 0.9572
bend 11 10 12 104.52
ENDCONSTRAINT
$end

$molecule
0 1
C -0.02800000086426735  0.004999999888241291  -0.0010000000474974513
H -0.4950000047683716  0.9340000152587891  0.3449999988079071
H -0.3440000116825104  -0.7979999780654907  0.671999990940094
H -0.3970000147819519  -0.210999995470047  -1.00600004196167
C 1.4859999418258667  0.1379999965429306  -0.020999999716877937
H 1.8669999837875366  0.34200000762939453  0.9900000095367432
H 1.9470000267028809  -0.7919999957084656  -0.3630000054836273
O 1.9299999475479126  1.1319999694824219  -0.9440000057220459
H 1.5729999542236328  1.9830000400543213  -0.6399999856948853
O -1.347000002861023  2.630000114440918  0.9760000109672546
H -0.9129999876022339  3.440000057220459  1.2410000562667847
H -2.2809998989105225  2.812000036239624  1.0820000171661377
$end

$coords
C -0.02800000086426735  0.004999999888241291  -0.0010000000474974513
H -0.4950000047683716  0.9340000152587891  0.3449999988079071
H -0.3440000116825104  -0.7979999780654907  0.671999990940094
H -0.3970000147819519  -0.210999995470047  -1.00600004196167
C 1.4859999418258667  0.1379999965429306  -0.020999999716877937
H 1.8669999837875366  0.34200000762939453  0.9900000095367432
H 1.9470000267028809  -0.7919999957084656  -0.3630000054836273
O 1.9299999475479126  1.1319999694824219  -0.9440000057220459
H 1.5729999542236328  1.9830000400543213  -0.6399999856948853
$end

This job does not converge even after 200 steps and the problem is that the Hessian stop updating after a certain point, with the following warning,

Cartesian Hessian Update
 Hessian updated using BFGS update
**WARNING** Hereditary positive definiteness endangered
 due to small product of step size times gradient change
**WARNING** Small scalar product
 Hessian update skipped this cycle

This is the gradient, displacement, and energy change at Step 200,

                            Maximum     Tolerance    Cnvgd?
         Gradient           0.006739      0.000300      NO
         Displacement       0.000002      0.001200     YES
         Energy change      0.000000      0.000001     YES

I have never encountered this problem before and not sure how to overcome this issue. I have tried increasing the tolerances (gradient up to 350 and displacement up to 1800) but it does not work.

#2

I can reproduce this behavior in Q-Chem v. 6.2.2. If the $opt section is removed (i.e., using the harmonic confining potentials but in an unconstrained optimization without the stretch and bend constraints), the optimization converges in 50 steps. Also curiously, using HF-3c (which I initially tried just for speed) the original, constrained optimization converges in 50 steps. Somehow the constraints are conspiring with the harmonic confiners to keep the gradient from going to zero but I don’t really understand it. Does the geometry that your are obtaining make sense to you?

#3

I also observed that the calculation converges by removing the water constraints, but that is not ideal since I need the water to stay as TIP3P water. The geometry at the end also does not make sense (compared to a similar job run in Orca), because the water oxygen should point towards the ethanol hydrogen, as shown in the image below (Orca output is yellow and Q-Chem is blue).

image
image794×730 61.2 KB