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 31 2019 - 23:20:45 CST