Geometry optimisation with constraints

Questions
#1

Dear Q-Chem users,

Is there a specific keyword that is required while performing constrained (spatially fixed/frozen atoms) optimization [DFT-based] in both the ground and the excited-state?
I end up exhausting the maximum optimization (energy seems to converge, I have given the energy of the last few optimization cycle below) cycle. I tried the following steps:

  1. Tweaking the geometries (changing the dihedrals of flexible portions of the molecule such as the aliphatic chain-substituents of aromatic ring) and reoptimize the structure.
  2. Reduce the default tolerance of the optimization parameters to a smaller value using “THRESH” keyword. I have also played around with specific keywords such as ( GEOM_OPT_TOL_GRADIENT , GEOM_OPT_TOL_ENERGY etc) but no success!
  3. I have also tried lengthening the geom_opt_cycles from default 50 to 400 cycles.

My rem-section had this
$rem
BASIS = cc-pVDZ
JOB_TYPE = Optimization
MEM_TOTAL = 120000
METHOD = wB97XD
SCF_CONVERGENCE = 6 ! reduced_convergence_criteria_for_thresh=9
SYMMETRY_IGNORE = 1
GEOM_OPT_MAX_CYCLES = 400
THRESH = 9 ! reduced_geom_convergence_criteria
$end

(Note: Thresh (and reduced scf convergence criteria) was used only after the default parameters didn’t give us ground state converged geometry)
$opt
FIXED
91 XYZ !HisC-spatially _fixed_atom
ENDFIXED
$end

Yet no success in the optimization of even the ground state of the molecule.

Energies of the last few optimization cycles:
-2707.349765079
-2707.349777164
-2707.349791392
-2707.349806339
-2707.349822085
-2707.349830906
-2707.349844057
-2707.349859100
-2707.349882210
-2707.349903446
-2707.349927979
-2707.349951354
-2707.349967897
-2707.349987905
-2707.350004919
-2707.350022825
-2707.350037582
-2707.350050833
-2707.350067663
-2707.350084984
-2707.350103316
-2707.350123358
-2707.350141690

Optimization parameters print at the last optimization cycle (400th) with Q-Chem reduced tolerance values for the optimization parameters:
Maximum Tolerance Cnvgd?
Gradient 0.001014 0.000300 NO
Displacement 0.053426 0.001200 NO
Energy change -0.000018 0.000001 NO

Its a fairly larger system with 112 atoms (tetrapyrrole ring with 2 neighboring 5-membered aromatic rings through hydrogen bonding interaction). The molecule is a closed-shell, +1 charged.

I was wondering if I am missing any constrained-optimization-specific keywords? I would also be happy to hear any suggestions you might have to get an optimized geometry with the constraints.

Thank you,

Siva

#2

Hello Siva,

So there are a couple things to note, first the THRESH keyword does not effect geometry optimization convergence criteria, it is associated with the cutoff for two electron integrals. This rem variable will only change the energy. The rem variables that will change convergence criteria are GEOM_OPT_TOL_GRADIENT, GEOM_OPT_TOL_DISPLACEMENT, and GEOM_OPT_TOL_ENERGY.

Also is there a specific reason you are trying to do constrained optimization by fixing a single atom? In general fixing a single atom is equivalent to doing an unconstrained optimization and then displacing the geometry to the fixed atom position. If for example you run water where the geometry convergence criteria and the SCF energy criteria are made extremely tight (for demonstration only, not in practice as these criteria greatly increase the job time)

$molecule
0 1
	O  0.000000000  0.000000000  0.369372944
	H -0.783975899  0.000000000 -0.184666472
	H  0.783975899  0.000000000 -0.184686472
$end

$rem
basis sto-3g
exchange hf
jobtype opt
scf_convergence 13
thresh 14
geom_opt_tol_gradient 1
geom_opt_tol_displacement 1
geom_opt_tol_energy 1
$end

Then re-run the job by fixing the oxygen

$opt
FIXED
1 XYZ
ENDFIXED
$end

You get the following results which differ by 2e-12 a.u.

Unconstrained: Final energy is   -74.9659011920794 

Constrained:   Final energy is   -74.9659011920773

A key reason I point this out is that in the geometry optimizer if fixed coordinates are specified the optimization can only be done in Cartesian space. Which in general is much slower to converge the geometry than using the default delocalized internal coordinates. This is also observed in the water example since the unconstrained optimization found the minimum in 6 cycles where the constrained optimization took 28 cycles to find the minimum.

I suggest that you remove the fixed atom position, rerun the geometry optimization with default settings thereby taking advantage of delocalized internal coordinate optimization. Then just manual displace your final structure to the original fixed atom position.

#3

Agree with Peter. A constrained optimization that fixes only one atom in space can pivot (rotate) around that one fixed atom in a manner that doesn’t change the energy (rigid rotation in free space) but is unlikely to ever satisfy the convergence criteria, except perhaps by accident.

#4

Thanks, @pfmclaug @jherbert.

@pfmclaug, I am trying to optimize just the QM region of the QM/MM setup. The QM region consists of a chromophore and few amino acids from the protein. To be able to relate the structures of the molecule and its protein neighbor in the QM region to its counterpart in the molecular dynamics trajectory, I need to freeze the frozen atoms. In the aforementioned example, I described only one frozen atom (which seemed more essential; to prevent the amino acid from moving around spatially a lot, which would not happen in a protein). I have also tried with frozen atoms on all the boundary atoms of the QM subsystem.

Thanks for the illustration with the water molecule, it was very useful. I did do unconstrained optimization (using delocalized coordinates) where I could get the equilibrium geometry but can’t relate with the structures in the MD trajectory (the amino acid that is otherwise frozen is way too flexible). The problem is, with 112 atoms the number of atoms in the QM system is already a lot (given we need also excited-state calculations to be performed) so adding few more amino acids in the QM region would be unfavorable concerning the computational costs.

I did manually displaced the equilibrium geometry of the unconstrained optimization but I get back to the same equilibrium geometry on reoptimization.

So I guess, using frozen atoms for optimizations may not be a great idea over the constraints with bonds/angles etc?

Is it possible to use (in optimizing an isolated QM subsystem i.e, without the environment) the scheme that is optimizing the QM subregion in the QM/MM optimization?

#5

For just this kind of purpose, my group developed harmonic confining potentials (https://doi.org/10.1021/acs.jpcb.9b11060). Take one of more atoms that you want to fix in place and place them in a harmonic well. It is described in the manual. On the other hand, the goal there was to be able to compute meaningful frequencies. If you don’t need to do that, then you might equally well be able to use fixed-atom constraints. My point was that a single fixed atom is going to lead to a situation where the molecule can just rotate around that fixed pivot point, changing the gradient but not the energy, in a manner that is unlikely to converge. I think this behavior should disappear as soon as you constrain at least two atoms.

#6

The problem is that a single atom being fixed leads to problems as demonstrated with water. Freezing multiple atoms will fix both the position and orientation of the molecule in space. Also freezing a bond is essentially equivalent to freezing two atoms. Therefore freezing atoms vs. freezing bonds/angles are in principal equivalent.

@jherbert suggestion of harmonic confined is a good one, but in general I think you need to fix additional atoms since you are under-defining the constraints.

#7

Thank you @jherbert. I will try and get back :slight_smile:

@pfmclaug, I also tried with multiple atoms frozen too but could not get it converged. As they were all boundary atoms of the QM region i.e, around 6 frozen atoms overall.

#8

We’ve had some similar issues with fixed-atom constrained geom opts that oscillate (when we cut active regions out of a protein crystal structure and tried to fix the C-alpha atoms to avoid collapse of the active site). This is the reason the atomic confining potentials were developed. It effectively accomplishes the same thing but in a “softer” way.

1 Like