Re: Bug in NAMD

From: Marcelo C. R. Melo (melomcr_at_gmail.com)
Date: Thu Jun 20 2019 - 09:46:36 CDT

Dear Pawel,

Thank you for the feedback and extensive tests.
In particular, your observation that the coordinates are repeated even with
"extraBonds off" and even using multiple QM packages.

Could you clarify:

1) In your different tests, were you always using a python wrapper script,
or did you ever try one of the native interfaces to ORCA or MOPAC?

2) How are you checking the coordinates sent to the QM software? Are you
saving a DCD at every step?

3) I noticed you are using the NAMD build from Nov/2018. (NAMD 2.13 for
Linux-x86_64-multicore, Built Fri Nov 9 14:39:16 CST 2018 by jim on
ganymede.ks.uiuc.edu)
Do you have access to a new one? Could you download the pre-compiled one
from the NAMD website for tests?
None of the issues solved in the newest releases would strike me as a
source for repeated coordinates being produced by the minimizer, but still,
other changes in NAMD could be affecting it.

Best,
Marcelo

On Thu, 20 Jun 2019 at 09:25, Pawel Kedzierski <pawel.kedzierski_at_pwr.edu.pl>
wrote:

> W dniu 18.06.2019 o 18:42, Brian Radak pisze:
>
> 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.
>
> No, GB only influences the MM part done by NAMD. I do have some water
> around the QM part but they also are not part of the QM subsystem, apart
> from being sent to QM program as point charges with electrostatic embedding.
>
> 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
>
> Yes, I am aware that the NAMD minimizers are only used to relax strains
> which would otherwise break dynamics. Still, I learned they are tuned
> against classical MM better than these higher order algorithms used by QM
> software and tuned for QM applications.
>
> The restraints, however, while important for my goals, are irrelevant for
> the problem I observe. With "extraBonds off" I still observe cycling of the
> same coordinates sent to the QM program with a stride of a few steps. And
> before I tested with constraints off, I used a range of force constants
> from 1000 down to 10.
>
> 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)
>
> True, but to some extent only. Well, the ultimate goal was to study
> dynamics effect and possibly also umbrella sampling.
>
> However being stuck with this NAMD problem I did switch to Gaussian for
> optimization, and so far I am rather disappointed. Neither the implicit
> solvent models nor the periodic boundary conditions in Gaussian are usable
> with molecular mechanics. Therefore, I solvated my system with a sphere of
> water and guess what? Optimization with Gaussian made the water explode
> away. Same system, same starting coordinates, same forcefield (Amber only,
> with QM part frozen) and in NAMD and water is optimized properly,
> preserving the spherical shape as expected due to surface tension.
>
> I guess that the Gaussian optimizer (which is second order) took steps too
> big what made the water clash first and then explode away, but I was
> unsuccessful with trials to limit the max step size taken and even when I
> started with water preoptimized with NAMD, in Gaussian it still would
> explode away.
>
> Also due to the performance of DFT/Amber QM/MM and frequent crashes
> ("error in internal coordinate system") I turned back to NAMD+ORCA.
>
> Thanks,
>
> Pawel
>
>
> 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 : Tue Dec 31 2019 - 23:20:49 CST