Questions about geometry optimization and frequency calculations with fixed atoms

Questions
#1

I was performing a geometry optimization + frequency calculation for acetate with the three H atoms fixed. Below is my input:

$molecule
-1 1
C          0.00000        0.00000        0.00000
C          0.00000        0.00000        1.55885
O         -1.04721        0.44384        2.08765
O          1.03473       -0.43856        2.11450
H         -0.40255       -0.94980       -0.36463
H         -0.66248        0.79074       -0.36465
H          1.00320        0.14979       -0.40377
$end

$rem
JOBTYPE  OPT 
METHOD  B3LYP
DFT_D  D3_BJ
BASIS  DEF2-TZVPD
SYM_IGNORE  TRUE
THRESH  14  
SCF_CONVERGENCE  8
MEM_TOTAL  8000
MEM_STATIC  2000
XC_GRID  000099000590
GEOM_OPT_TOL_GRADIENT  100 
GEOM_OPT_TOL_DISPLACEMENT  100 
GEOM_OPT_TOL_ENERGY  10  
$end

$opt
 FIXED
 5:7  XYZ 
 ENDFIXED
$end

@@@

$molecule
read
$end

$rem
JOBTYPE  FREQ
METHOD  B3LYP
DFT_D  D3_BJ
BASIS  DEF2-TZVPD
SYM_IGNORE  TRUE
THRESH  14  
SCF_CONVERGENCE  8
MEM_TOTAL  8000
MEM_STATIC  2000
XC_GRID  000099000590
FRZN_OPT  TRUE
FRZ_ATOMS 3
$end

$frozen_opt
5 6 7
$end

What I observed:

  1. The positions of the three H atoms are changing during the optimization
  2. With the “zeroing out” Hessian approach applied, the normal modes still involve contributions from the three H atoms.
 Mode:                 1                      2                      3
 Frequency:       800.95                1316.43                1387.54
 Force Cnst:      4.0235                 4.7675                 4.2042
 Red. Mass:      10.6451                 4.6692                 3.7064
 IR Active:          YES                    YES                    YES
 IR Intens:       62.028                165.375                104.249
 Raman Active:       YES                    YES                    YES
               X      Y      Z        X      Y      Z        X      Y      Z    
 C          0.001  0.005 -0.037    0.006  0.010  0.294   -0.010 -0.025 -0.317
 C         -0.036 -0.056 -0.481    0.001 -0.021  0.249    0.011  0.007  0.329
 O          0.432 -0.156  0.155    0.181 -0.064 -0.178    0.099 -0.041 -0.037
 O         -0.412  0.181  0.179   -0.183  0.081 -0.164   -0.104  0.040 -0.035
 H          0.050  0.099  0.374   -0.151 -0.370 -0.351    0.180  0.452  0.303
 H          0.104  0.018  0.241   -0.216  0.215 -0.334    0.217 -0.206  0.338
 H         -0.058  0.077  0.251    0.321  0.015 -0.354   -0.324 -0.007  0.356
 TransDip  -0.007 -0.000 -0.252    0.007 -0.008  0.412    0.007 -0.017  0.326

 Mode:                 4                      5                      6
 Frequency:      1570.71                1635.91                1709.55
 Force Cnst:      2.0656                19.1554                 3.1420
 Red. Mass:       1.4210                12.1485                 1.8247
 IR Active:          YES                    YES                    YES
 IR Intens:       27.226                702.139                 54.665
 Raman Active:       YES                    YES                    YES
               X      Y      Z        X      Y      Z        X      Y      Z    
 C         -0.063 -0.162  0.029   -0.048  0.019  0.001    0.226 -0.089  0.001
 C          0.028  0.071 -0.012    0.786 -0.313 -0.024   -0.035  0.014 -0.001
 O         -0.007 -0.011 -0.004   -0.279  0.106  0.140   -0.027  0.010  0.065
 O         -0.003 -0.013 -0.004   -0.281  0.118 -0.123   -0.029  0.013 -0.065
 H          0.264  0.656  0.004    0.007  0.000  0.001   -0.220  0.073  0.000
 H         -0.033  0.477 -0.038    0.032 -0.036 -0.174   -0.513  0.428 -0.084
 H          0.348  0.330 -0.044    0.052  0.002  0.173   -0.642  0.033  0.097
 TransDip  -0.061 -0.155  0.018    0.788 -0.315 -0.023    0.220 -0.087 -0.005

I am thus not sure if this calculation is working correctly.

#2

Hi Yuezhi, I checked the output of the simpler Manual Example 9.12 with v6.2.1 (9.4.3 Frozen Atoms‣ 9.4 Constrained Optimization ‣ Chapter 9 Exploring Potential Energy Surfaces: Searches for Critical Points and Molecular Dynamics ‣ Q-Chem 6.2 User’s Manual) and although the coordinates of the fixed atoms change during the optimization, their relative positions don’t. This is also true in your OPT calculation. I have not yet looked at the hessian contribution though.

#3

Yes, you are right. I also checked the relative positions between the three H atoms in my calculation, and noticed the distances between them remain the same during the optimization. So I guess the fixed-atom optimization might be fine.

The resulting normal modes from the frequency job is more puzzling to me right now.

#4

I am not sure what exactly you mean by the “zeroing out hessian” approach, so please correct me. In the context of optimization, say with methods that involve an approximate hessian such as BFGS, projecting out the contribution to the gradient from these hydrogens would be equivalent to zeroing out the rows and columns of the approximate hessian corresponding to the coordinates of these hydrogens. But, that has got nothing to do with the frequency calculation, correct? The electronic hessian in the context of a frequency calculation and its eigenvectors should have contributions from the hydrogens, fixed or not.

#5

The keywords FRZN_OPT = TRUE and FRZ_ATOMS = 3 (plus the $frozen_opt section) in the 2nd job actually has an effect on the vibrational frequency analysis, which is documented on the manual: 9.4.3 Frozen Atoms‣ 9.4 Constrained Optimization ‣ Chapter 9 Exploring Potential Energy Surfaces: Searches for Critical Points and Molecular Dynamics ‣ Q-Chem 6.2 User’s Manual