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

From: Brian Radak (brian.radak_at_gmail.com)
Date: Tue Jun 18 2019 - 11:42:02 CDT

Does NAMD QM work with GB? I'm not aware of any mechanism by which that
would enter the QM calculation in the packages you listed.

If you are just doing optimization:

1) for a large system you are quite unlikely to find any true minima - the
PESs are just too noisy. Your restraints may also still be too strong

2) unless you plan to perform dynamics next, there probably isn't much
utility in using NAMD for this, since the QM packages should be able to
handle this natively (and with much better minimizers)

BKR

On Tue, Jun 18, 2019, 2:49 AM Pawel Kedzierski <pawel.kedzierski_at_pwr.edu.pl>
wrote:

> 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 "
> qmForces on".
>
> 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 : Wed Sep 18 2019 - 23:21:03 CDT