Big oscillations during QM/MM minimization

From: Pawel Kedzierski (pawel.kedzierski_at_pwr.edu.pl)
Date: Tue May 07 2019 - 08:53:18 CDT

Dear NAMD experts,

I am trying to employ namd for QMMM optimization with ORCA to model an
enzymatic reaction. To mimic a relaxed scan along the reaction
coordinate, I use harmonic constraints (ExtraBonds) to enforce a
specific length of a distance in the QM region.

During energy minimization, I observe regular behavior of the minimizer
for the first few hundred steps, but then I get oscillations of the
energy values of constant amplitude and precisely periodic. The
amplitude is of the order of 100 kcal/mole, and the energies are
repeated verbatim every few steps:

QMENERGY:    1492         1.0000  -1248781.5890 <-
QMENERGY:    1493         1.0000  -1248664.2252
QMENERGY:    1494         1.0000  -1248729.0129
QMENERGY:    1495         1.0000  -1248728.6078
QMENERGY:    1496         1.0000  -1248781.5890 <-
QMENERGY:    1497         1.0000  -1248664.2252
QMENERGY:    1498         1.0000  -1248729.0129
QMENERGY:    1499         1.0000  -1248728.6078
QMENERGY:    1500         1.0000  -1248781.5890 <-

And the total energies:

awk '/^ENERGY:/{print $2, $12}' QMMM_2.0.log | tail -n 9
1492 -1267497.5702 <-
1493 -1267319.3505
1494 -1267444.5801
1495 -1267444.6028
1496 -1267497.5702 <-
1497 -1267319.3505
1498 -1267444.5801
1499 -1267444.6028
1500 -1267497.5702 <-

Note that when I only change "qmForces" to "off", the oscillations do
not happen and minimization goes on normally.

What I have tried, without effect so far:

  * lowering the force constant of the ExtraBonds constraint: 1000, 500, 200
  * lowering minBabyStep to 0.0001
  * changing the minimizer to velocityQuenching with maximumMove 0.0001

None of these helped. What else can I try?

I also saved the forces using "output onlyforces file.pdb" and the force
components are nowhere near zero but rather in the range 1-50 for most
of the atoms (bot QM and MM parts). This is another surprise for me. I
would expect the forces to approach zero. Even if there is a problem
with QM part, most of the MM part should be relaxed, and for the
crystallographic water molecules I think this should be rather quick.

Am I missing something here?

Constrained distance file:

# Reaction coordinate: O - P distance
bond 7324 7289 1000.0 2.0

Config excerpt (using CHARMM forcefield version 36):

cutoff 16
pairlistdist 18.0
switching on
switchdist 15
exclude scaled1-4
1-4scaling 1.0
rigidbonds none

solventDielectric 80.0
sasa on
gbis on
alphaCutoff 14.0
ionConcentration 0.15

QMBondScheme "CS"
QMSwitching on
QMSwitchingType Switch
QMPointChargeScheme Round
qmReplaceAll OFF
QMVdWParams off
qmElecEmbed On
QMPCStride 1
QMCustomPCSelection off
QMLiveSolventSel off
# ORCA options
qmConfigLine  "! HF-3c ENGRAD"
qmConfigLine  "%%output PrintLevel Mini Print\[P_Mulliken\] 1 Print\[P_AtCharges
qmConfigLine  "%%pal nprocs 4 end"

Versions:

NAMD 2.13 for Linux-x86_64-multicore, Built Fri Nov 9 14:39:16 CST 2018
by jim on ganymede.ks.uiuc.edu

ORCA static binary release version 4.0.1.2 with OpenMPI 2.1.2

Gentoo GNU/Linux 4.1.12 on Intel Xeon workstation.

With regards,

Paweł Kędzierski


This archive was generated by hypermail 2.1.6 : Tue Dec 10 2019 - 23:20:45 CST