Bug in NAMD (was: Big oscillations during QM/MM minimization)

From: Pawel Kedzierski (pawel.kedzierski_at_pwr.edu.pl)
Date: Tue Jun 18 2019 - 01:47:24 CDT

Dear NAMD experts,

I have adapted the script runORCA.py distributed with
NAMD_2.13_Linux-x86_64 in the subdirectory lib/qmmm to see, what is
going on during QM/MM minimization.
The observation is, that *every few steps NAMD saves the exactly same
coordinates for the QM program to calculate forces* (my script compares
them with epsilon 2e-6 as the precision in the file written by NAMD is
to the sixth decimal point). This behavior is consistent independently
from the QM program used (MOPAC, ORCA or CUSTOM), the minimizer (cg or
velocity quenching) or any other QM options I tried to set in the config
file. Even if I modify the forces returned from the QM calculations,
still the program keeps to cycle the same set of coordinates again and
again.
Identical coordinates of the QM are not saved every time, but this
behavior start to appear early on:
REPEATED COORDS from step000003 at step 6

and later i see:
REPEATED COORDS from step000003 at step 13
REPEATED COORDS from step000006 at step 13
REPEATED COORDS from step000009 at step 13

and at the end:
.. cut out ..
REPEATED COORDS from step004977 at step 5001
REPEATED COORDS from step004980 at step 5001
REPEATED COORDS from step004983 at step 5001
REPEATED COORDS from step004986 at step 5001
REPEATED COORDS from step004989 at step 5001
REPEATED COORDS from step004992 at step 5001
REPEATED COORDS from step004995 at step 5001
REPEATED COORDS from step004998 at step 5001

The stride is also not always 3. There is something very wrong, no
minimization algorithm should behave like this; and it only happens with
"qmForceson".

With regards,
Paweł Kędzierski

W dniu 07.05.2019 o 15:53, Pawel Kedzierski pisze:
>
> 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 : Thu Dec 05 2019 - 23:20:45 CST