From: Pawel Kedzierski (pawel.kedzierski_at_pwr.edu.pl)
Date: Thu Jun 20 2019 - 10:35:11 CDT
W dniu 20.06.2019 o 16:46, Marcelo C. R. Melo pisze:
> 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?
No, I started using the native interfaces (both of them) and switched to
the runORCA.py script to actually identify the source of the problem.
>
> 2) How are you checking the coordinates sent to the QM software? Are
> you saving a DCD at every step?
Well, it did not occur to me until now to compare the DCD frames, but I
guess they should show the same behaviour.
I used numpy in runORCA.py to convert the coordinates received from NAMD
to numpy array, and save these coordinates in a file.
Every time the script is executed by NAMD it reads the tables of all
coordinates saved previously, keyed with step number, and compares them
to the current ones:
import numpy as np
from os.path import exists
prevEps = 2e-6
if exists(prevFile):
tmp = np.load(prevFile)
prevCoord = dict(tmp)
tmp.close()
steps = sorted([k for k in prevCoord.keys() if k[:4] == 'step'])
last = int(steps[-1][4:])
for s in steps:
if np.all(np.abs(QMcoords - prevCoord[s]) < prevEps):
print("REPEATED COORDS from %s at step %d" % (s, last+1))
....
The script is in available within this zipped example:
http://mml.ch.pwr.wroc.pl/p.kedzierski/namd/implicit_CHARMM.zip
>
> 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 <http://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?
Yes, I could try this out.
With regards,
Paweł
> 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 <mailto: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
>> <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 Dec 31 2020 - 23:17:11 CST