I am trying to do a constrained optimization with constraint on an out-of-plane angle. Input structure satisfies the constraint. As the optimization cycles proceed, the angle constraint is increasingly violated. This could be a mistake on my part, but the frozen_bend and frozen_dihedral work fine for me with the same syntax.
Background: We want to study the interaction energy between a dimer consisting of pyridine and N-methyl acetamide. This dimer forms a hydrogen bond between the N(pyr) and the H(nma)-N(nma) and we want to see how the nature of the interaction changes as the Hydrogen on the NMA is taken out of the plane of the pyridine ring. The plane of the pyridine ring is defined by the Nitrogen and its two adjacent Carbon sites. The angle we want to constrain is the angle between this ring and a vector that goes from the N(pyr) to the hydrogen bonding H(nma).
Here is an input file exhibiting the problem:
import psi4
psi4.set_memory(‘10 GB’)
set G_convergence gau_tight
set geom_maxiter 500
set geom_maxiter 5
set maxiter 300
set scf_type DF
set freeze_core true
set e_convergence 11
set d_convergence 11
set r_convergence 11
set ints_tolerance 0.0
set df_basis_scf aug-cc-pvtz-jkfit
set df_basis_mp2 aug-cc-pvtz-ri
pyrid_nma3 = psi4.geometry(“”"
#xi 60.00 ; nu 20.00 ; s 1.00 ; R( 6-19) 1.9879 ; phi(19- 6- 5- 7) 109.88 ; theta(19- 6- 5) 112.94
0 1
H 4.178003 1.272543 -1.390568
H 1.800200 1.236473 -2.194888
C 3.676151 0.359073 0.526447
C 3.388898 0.867726 -0.751420
C 2.058413 0.848600 -1.204580
N 1.031900 0.362451 -0.470541
C 1.319026 -0.125488 0.759040
C 2.616476 -0.147780 1.297562
H 0.466918 -0.513626 1.327423
H 2.786709 -0.554941 2.297623
H 4.699761 0.357627 0.911431
0 1
C -2.425910 0.569105 0.634079
H -3.367442 0.700810 0.081109
H -1.670000 0.127229 -0.031852
H -2.624419 -0.112909 1.474234
C -1.999320 1.918163 1.192212
O -2.713116 2.564496 1.975950
N -0.784426 2.369788 0.770789
H -0.204988 1.781993 0.167202
C -0.269249 3.641140 1.257976
H -0.172366 3.633280 2.355169
H -0.946572 4.464895 0.986561
H 0.716058 3.813399 0.803668
units angstrom
symmetry c1
noreorient
“”")
set opt_coordinates cartesian
custom_opt_params = {
“H-bond_connect”: True ,
“frozen_oofp”: “19 6 7 5”,
}
energy = psi4.optimize(‘MP2/aug-cc-pvdz’, molecule=pyrid_nma3,
engine=‘optking’,
optimizer_keywords=custom_opt_params)
print(f"Final MP2 energy: {energy} Hartree")
Here is some of the out-of-plane trajectory from the output file. It continues to drift:
*O(19,6,7,5) = 1.047198 60.000054
*O(19,6,7,5) = 1.055156 60.455985
*O(19,6,7,5) = 1.053929 60.385676
*O(19,6,7,5) = 1.049238 60.116936
*O(19,6,7,5) = 1.049878 60.153556
*O(19,6,7,5) = 1.058147 60.627341
*O(19,6,7,5) = 1.052421 60.299255
*O(19,6,7,5) = 1.048358 60.066499
*O(19,6,7,5) = 1.049118 60.110021
*O(19,6,7,5) = 1.047951 60.043164
*O(19,6,7,5) = 1.043915 59.811938
*O(19,6,7,5) = 1.030272 59.030260
*O(19,6,7,5) = 1.031478 59.099350
*O(19,6,7,5) = 1.025006 58.728503
*O(19,6,7,5) = 1.024342 58.690449
*O(19,6,7,5) = 1.021760 58.542515
*O(19,6,7,5) = 1.020364 58.462523
*O(19,6,7,5) = 1.015939 58.209043
*O(19,6,7,5) = 1.011254 57.940570
*O(19,6,7,5) = 1.001806 57.399279
*O(19,6,7,5) = 1.000464 57.322345