QM/MM energy conservation
139 views
Skip to first unread message
mayank...@gmail.com
Aug 4, 2020, 2:52:01 AM8/4/20
to cp2k
Hi,
I am trying a NVT MD simulation for a QM/MM system with QM water molecules surrounding a spatially fixed MM molecule, which are in turn surrounded by MM water molecules. I am using a time step of 5 a.u. (for particular reasons) and the QM box is much smaller than compared to full MM box so within the timescale of MD simulation the QM molecules do not cross the QM box boundary. I am attaching the relevant sections of input file:
----
&MOTION
&MD
ENSEMBLE NVT
STEPS 10000
TIMESTEP 1.21E-01
TEMPERATURE 3.00E+02
&THERMOSTAT
TYPE NOSE
REGION DEFINED
&DEFINE_REGION
QM_SUBSYS ATOMIC
&END DEFINE_REGION
&DEFINE_REGION
MM_SUBSYS ATOMIC
&END DEFINE_REGION
&NOSE
TIMECON 9.99E+02
&END NOSE
&END THERMOSTAT
&FORCE_EVAL
METHOD QMMM
&DFT
BASIS_SET_FILE_NAME ./GTH_BASIS_SETS
POTENTIAL_FILE_NAME ./GTH_POTENTIALS
CHARGE 0
&SCF
MAX_SCF 60
EPS_SCF 2.0E-07
SCF_GUESS RESTART
&OT T
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&OUTER_SCF T
EPS_SCF 2.0E-07
MAX_SCF 40
&END OUTER_SCF
&PRINT
&RESTART SILENT
BACKUP_COPIES 1
&END RESTART
&END PRINT
&END SCF
&QS
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 4
&END QS
&MGRID
CUTOFF 3.2E+02
COMMENSURATE T
&END MGRID
&XC
DENSITY_CUTOFF 1.0E-10
GRADIENT_CUTOFF 1.0E-10
TAU_CUTOFF 1.0E-10
&XC_GRID
XC_SMOOTH_RHO NN50
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&XC_FUNCTIONAL NO_SHORTCUT
&PBE T
PARAMETRIZATION REVPBE
&END PBE
&END XC_FUNCTIONAL
&VDW_POTENTIAL
POTENTIAL_TYPE PAIR_POTENTIAL
&PAIR_POTENTIAL
R_CUTOFF 1.9E+01
TYPE DFTD3
PARAMETER_FILE_NAME ./dftd3.dat
REFERENCE_FUNCTIONAL revPBE
CALCULATE_C9_TERM T
REFERENCE_C9_TERM T
LONG_RANGE_CORRECTION F
&END PAIR_POTENTIAL
&END VDW_POTENTIAL
&END XC
&END DFT
&MM
&FORCEFIELD
...
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE SPME
RCUT 1.2E+01
ALPHA 2.917E-01
GMAX 75
O_SPLINE 4
&END EWALD
&END POISSON
&END MM
&QMMM
E_COUPL GAUSS
USE_GEEP_LIB 6
NOCOMPATIBILITY T
CENTER NEVER
INITIAL_TRANSLATION_VECTOR 0.0E+00 0.0E+00 0.0E+00
&QM_KIND O
MM_INDEX ...
&END QM_KIND
&QM_KIND H
MM_INDEX ...
&END QM_KIND
&MM_KIND H
RADIUS 4.4E-01
&END MM_KIND
&MM_KIND O
RADIUS 1.2E+00
&END MM_KIND
.....
&CELL
ABC 4.0E+01 3.0E+01 3.0E+01
PERIODIC XYZ
&END CELL
&PERIODIC
GMAX 5.0E-01
&MULTIPOLE ON
RCUT 1.2E+01
EWALD_PRECISION 1.0E-08
ANALYTICAL_GTERM T
&END MULTIPOLE
&END PERIODIC
&END QMMM
&SUBSYS
&CELL
A 5.0E+01 0.0E+00 0.0E+00
B 0.0E+00 4.0E+01 0.0E+00
C 0.0E+00 0.0E+00 4.0E+01
MULTIPLE_UNIT_CELL 1 1 1
&END CELL
&COORD
-------
With this setup for the production run I obtained the following trends for the potential and constant quantity as the graphs attached here. I have tested the QM setup for a fully QM calculation on bulk water, which works fine. My main aim here is to reduce the energy drift for the QM/MM system, so is there any other aspect of the calculation should I take care to reduce the drift? Any suggested corrections for the QM/MM file would also be helpful.
Best Regards,
Mayank Dodia
Thomas Kühne
Aug 5, 2020, 4:57:17 AM8/5/20
to 'Dorothea Golze' via cp2k
Dear Mayank,
the y-axis of your conserved quantity plot is too coarse to judge if energy conservation is an issue at all.
Possibly it is already good enough, but in any case the integration time step error is not s.th. you should
be concerned first and foremost, hence you can safely increase TIMESTEP. Also, since energy is an extensive
property, in my opinion the best way to look at the conserved quantity is in units of Kelvin per atom and ps.
The potential energy is a matter of your initial configuration and only shows that you haven’t equilibrated
as yet …
Cheers,
Thomas
--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/63b0f162-c684-481d-884a-f10890c1aa91n%40googlegroups.com.
<const_quantity.png><pot_enr.png>
==============================
Thomas D. Kühne
Dynamics of Condensed Matter
Chair of Theoretical Chemistry
University of Paderborn
Warburger Str. 100
D-33098 Paderborn
Germany
Mayank Dodia
Aug 5, 2020, 5:19:02 AM8/5/20
to cp2k
Dear Thomas,
Thanks for your response. My aim is to understand the solvation shell around the central MM molecule, so I kept the timestep and the equilibration length short to avoid excessive motion of QM water molecules and decomposing the solvation shell. Do you have any recommendations on how to accelerate equilibration in such a scenario?
Best,
Mayank
On Wednesday, August 5, 2020 at 10:57:17 AM UTC+2, tkuehne wrote:
Dear Mayank,the y-axis of your conserved quantity plot is too coarse to judge if energy conservation is an issue at all.Possibly it is already good enough, but in any case the integration time step error is not s.th. you shouldbe concerned first and foremost, hence you can safely increase TIMESTEP. Also, since energy is an extensiveproperty, in my opinion the best way to look at the conserved quantity is in units of Kelvin per atom and ps.The potential energy is a matter of your initial configuration and only shows that you haven’t equilibratedas yet …Cheers,Thomas
To unsubscribe from this group and stop receiving emails from it, send an email to cp...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/63b0f162-c684-481d-884a-f10890c1aa91n%40googlegroups.com.
<const_quantity.png><pot_enr.png>
Reply all
Reply to author
Forward
0 new messages