Re: Bug in NAMD

From: Pawel Kedzierski (pawel.kedzierski_at_pwr.edu.pl)
Date: Thu Jun 20 2019 - 08:23:43 CDT

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 <mailto: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 "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 <http://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 Sep 19 2019 - 23:20:36 CDT