Dear AMBER Developers,
I would like to report a potential bug in the virial calculation during QM/MM simulations in AMBER.
In my test, I simulated a system of 489 methanol molecules using a classical MM force field. I then converted one methanol molecule to the QM region using DFTB3 (via DFTB+), keeping the rest as MM. I enabled verbose=3 in the Ewald namelist to monitor the virial components.
Upon comparison, I observed that most virial tensor terms changed only slightly between the MM and QM/MM simulations. However, the Reciprocal VIRIAL (rec_vir in ew_force.F90) showed a substantial discrepancy: in the pure MM simulation, the diagonal elements of rec_vir are on the order of ~1e3, while in the QM/MM case, they are nearly zero. This discrepancy appears to propagate through the Atomic VIRIAL and eventually affects the Molecular VIRIAL, which is used to compute the system pressure.
>From my reading of the code (ew_force.F90), I suspect that when qm_pme is enabled, rec_vir is incorrectly computed using only the QM charges. This could potentially result in a lower equilibrium density in QM/MM simulations.
Due to the complexity of the code, I haven’t been able to trace the issue fully. I’ve attached my output files (md_qmmm.out) from both the MM-only and QM/MM simulations (in folders MOH_489mm_0qm and MOH_488mm_1qm) for your reference. I’ve also included excerpts of the relevant output. If internal developer documentation exists, I would appreciate access, as it would greatly assist my understanding and debugging. I’d also be interested in joining the developer list to contribute directly, if possible.
Additionally, I’d like to make a suggestion regarding the Ewald error estimate (virvsene in ew_force.F90). In QM/MM simulations, the estimate is currently around 0.1 ( in contrast to values near 1e-5 seen in pure MM simulations ), which I understand may not be a meaningful error because it includes inconsistent contributions (QM-QM, QM-MM, MM-MM). It would be helpful if this estimate could be modified to include only the QM-MM and MM-MM parts to ensure consistency. This would make the error estimate more useful for convergence testing (of QM/MM Ewald/PME parameters) and debugging in QM/MM runs.
Thank you for your time and for developing such a powerful simulation package.
Best regards,
Zihao Bai
link to the attached files:
[
https://res.public.onecdn.static.microsoft/assets/mail/file-icon/png/generic_16x16.png]amber_bug.tar.gz<
https://uwprod-my.sharepoint.com/:u:/g/personal/zbai29_wisc_edu/EfOqce-qoD1IlVl8DjxQcw8BHNsGmKCY1eKF23TLyUYwZg>
MOH_488mm_1qm/md_qmmm.out
MOLECULAR VIRIAL: -.19264834E+03 0.71709577E+02 -.42368592E+02
MOLECULAR VIRIAL: -.67195242E+01 -.18109503E+03 0.20552265E+03
MOLECULAR VIRIAL: 0.94223511E+01 0.12618927E+03 -.67950775E+02
--------------------------------------------
Reciprocal VIRIAL: 0.59187371E+00 -.59514589E-01 -.31221085E+00
Reciprocal VIRIAL: -.59514589E-01 0.16435798E+00 -.44681210E+00
Reciprocal VIRIAL: -.31221085E+00 -.44681210E+00 0.15767102E+01
..................
Direct VIRIAL: -.18771602E+04 -.11018011E+02 0.50294480E+02
Direct VIRIAL: -.11018011E+02 -.19356208E+04 0.15414421E+03
Direct VIRIAL: 0.50294480E+02 0.15414421E+03 -.21494256E+04
..................
Dir Sum EE vir trace: 0.39380183E+04
..................
Adjust VIRIAL: -.41850648E+03 -.40964540E+01 -.19596137E+01
Adjust VIRIAL: -.40964540E+01 -.40392793E+03 0.12612842E+02
Adjust VIRIAL: -.19596137E+01 0.12612842E+02 -.42534860E+03
..................
Recip Disp. VIRIAL: 0.25181754E+03 0.00000000E+00 0.00000000E+00
Recip Disp. VIRIAL: 0.00000000E+00 0.25181754E+03 0.00000000E+00
Recip Disp. VIRIAL: 0.00000000E+00 0.00000000E+00 0.25181754E+03
..................
Self VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Self VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Self VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
..................
E14 VIRIAL: -.88327310E+03 0.86200384E+01 0.13350157E+02
E14 VIRIAL: 0.86200384E+01 -.86675092E+03 0.28609997E+02
E14 VIRIAL: 0.13350157E+02 0.28609997E+02 -.84812120E+03
..................
Atomic VIRIAL: -.29265304E+04 -.65539410E+01 0.61372812E+02
Atomic VIRIAL: -.65539410E+01 -.29543178E+04 0.19492024E+03
Atomic VIRIAL: 0.61372812E+02 0.19492024E+03 -.31695012E+04
--------------------------------------------
Sub VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Sub VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Sub VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
--------------------------------------------
NSTEP = 0 TIME(PS) = 1005.000 TEMP(K) = 305.87 PRESS = 908.0
Etot = -2282.9100 EKtot = 2080.5896 EPtot = -4363.4997
BOND = 151.6774 ANGLE = 909.1122 DIHED = 141.8674
1-4 NB = 0.0000 1-4 EEL = 2598.1452 VDWAALS = -908.0639
EELEC = -3152.1618 EHBOND = 0.0000 RESTRAINT = 0.0000
DFTBPLUSESCF= -4104.0762
EKCMT = 410.9522 VIRIAL = -220.8471 VOLUME = 32227.7558
Density = 0.8069
Ewald error estimate: 0.1458E+00
------------------------------------------------------------------------------
MOH_489mm_0qm/md_qmmm.out
MOLECULAR VIRIAL: -.54250569E+02 0.63261131E+02 -.33843649E+02
MOLECULAR VIRIAL: -.16352381E+02 -.17710853E+02 0.18728911E+03
MOLECULAR VIRIAL: 0.18970749E+02 0.10769364E+03 0.82787875E+02
--------------------------------------------
Reciprocal VIRIAL: 0.15099890E+03 -.93260004E+01 0.75365558E+01
Reciprocal VIRIAL: -.93260004E+01 0.16473493E+03 -.10266706E+02
Reciprocal VIRIAL: 0.75365558E+01 -.10266706E+02 0.14578722E+03
..................
Direct VIRIAL: -.18764924E+04 -.12328384E+02 0.50732170E+02
Direct VIRIAL: -.12328384E+02 -.19285661E+04 0.15458958E+03
Direct VIRIAL: 0.50732170E+02 0.15458958E+03 -.21373047E+04
..................
Dir Sum EE vir trace: 0.39578618E+04
..................
Adjust VIRIAL: -.41913593E+03 -.40416912E+01 -.16378436E+01
Adjust VIRIAL: -.40416912E+01 -.40412820E+03 0.13156109E+02
Adjust VIRIAL: -.16378436E+01 0.13156109E+02 -.42713890E+03
..................
Recip Disp. VIRIAL: 0.25181754E+03 0.00000000E+00 0.00000000E+00
Recip Disp. VIRIAL: 0.00000000E+00 0.25181754E+03 0.00000000E+00
Recip Disp. VIRIAL: 0.00000000E+00 0.00000000E+00 0.25181754E+03
..................
Self VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Self VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Self VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
..................
E14 VIRIAL: -.88746220E+03 0.96747666E+01 0.13199460E+02
E14 VIRIAL: 0.96747666E+01 -.86759156E+03 0.28559233E+02
E14 VIRIAL: 0.13199460E+02 0.28559233E+02 -.84852098E+03
..................
Atomic VIRIAL: -.27802741E+04 -.16021309E+02 0.69830342E+02
Atomic VIRIAL: -.16021309E+02 -.27837334E+04 0.18603822E+03
Atomic VIRIAL: 0.69830342E+02 0.18603822E+03 -.30153598E+04
--------------------------------------------
Sub VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Sub VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
Sub VIRIAL: 0.00000000E+00 0.00000000E+00 0.00000000E+00
--------------------------------------------
NSTEP = 0 TIME(PS) = 1005.000 TEMP(K) = 305.87 PRESS = 582.8
Etot = 1810.9117 EKtot = 2080.5896 EPtot = -269.6779
BOND = 151.7204 ANGLE = 909.9810 DIHED = 141.9360
1-4 NB = 0.0000 1-4 EEL = 2603.5747 VDWAALS = -908.0639
EELEC = -3168.8261 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 410.9522 VIRIAL = 5.4132 VOLUME = 32227.7558
Density = 0.8069
Ewald error estimate: 0.4850E-04
------------------------------------------------------------------------------
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Jun 28 2025 - 09:30:02 PDT
