Re: Big oscillations during QM/MM minimization

From: Dawid das (addiw7_at_googlemail.com)
Date: Tue May 07 2019 - 10:07:20 CDT

Hi Paweł,

I have never performed QM/MM calculations with ORCA/NAMD but as far as I
understand you use electrostatic embedding.
Obviously, it is important, but did you give a try with mechanical
embedding? If the oscillations won't appear it will perhaps
give you an idea about what is going on there.

Best wishes,
Dawid Grabarek

wt., 7 maj 2019 o 15:54 Pawel Kedzierski <pawel.kedzierski_at_pwr.edu.pl>
napisał(a):

> 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 : Fri Dec 06 2019 - 23:20:38 CST