Cons Qty[a.u.] drift? kinetic + potential = Cons Qty ?
marco.stenta
marco.... at gmail.com
Thu Aug 4 07:06:38 UTC 2011
Dear all,
I have a question for any QM/MM user/developer.
I have a NVT simulation of a big protein embedded in a water box,
already equilibrated at MM level, and then 2 ps were performed by
keepign the heavy QM atoms frozen, as to further relax the MM atoms
surrounding the active site.
*the system contains a couple of Mg atoms (Mg: TZVDD3DF3PD-GTH-BLYP
all the other atoms DZVP-GTH-BLYP)
*charge of the qm system is -1, the spurious charge left on the MM
part was somehow redistributed/adjusted as to have
CHARGE_INFO| Total Charge of the Classical
System: 0.000006
*QM is not periodic and MT is used as poisson solver for the QM part
(the box was chosen to be 2.5 times the size of the QM part, which is
globular and compact)
* the cutoff was set to 320 and the grid is COMMENSURATE
*the parameter XC_SMOOTH_RHO NN50 was also set
*H are substituted with D, and the timestep is set to 0.2 fs
(the full input is attached below)
what troubles me is the drift in Cons Qty.
I understand I can still be out of equilibrium and the system should
further relax, but in similar simulations I observed a constant drift.
Moreover the Cons Qty[a.u.] is not the sum of potential and kinetic
energy (.ener file): is this OK? or it is a bell ringing for
something deeply wrong in my calculation/setup/input/system?
cons qty -(pot + kin) = err
step 1: -1517.083028953-(-1747.938628518+230.841109714) = .
014489851
step 203: -1517.034195266-(-1747.682577807+230.042654659) = .
605727882
Thansk a lot
Marco
# Step Nr. Time[fs] Kin.[a.u.]
Temp[K] Pot.[a.u.] Cons Qty[a.u.] UsedTime[s]
0 0.000000 230.867136550
300.555492024 -1747.950286957 -1517.083150407
0.000000000
1 0.200000 230.841109714
300.521608861 -1747.938628518 -1517.083028953
399.900000000
2 0.400000 230.807655637
300.478056511 -1747.928456404 -1517.083062027
134.000000000
3 0.600000 230.730159995
300.377168437 -1747.920690429 -1517.082994355
130.840000000
4 0.800000 230.740950337
300.391215894 -1747.914613364 -1517.081103127
130.830000000
5 1.000000 230.761638258
300.418148564 -1747.914568012 -1517.081032391
134.510000000
6 1.200000 230.746886303
300.398943661 -1747.920038506 -1517.081850762
134.440000000
7 1.400000 230.781912025
300.444542065 -1747.929605498 -1517.082015560
130.700000000
.
.
.
.
.
198 39.600000 230.304901871
299.823544102 -1747.846813641 -1517.033160667
132.440000000
199 39.800000 230.235992906
299.733834635 -1747.807562871 -1517.033321937
122.030000000
200 40.000000 230.132127496
299.598616952 -1747.770362015 -1517.033673946
125.460000000
201 40.200000 230.145587240
299.616139584 -1747.737129241 -1517.034724628
134.050000000
202 40.400000 230.026559444
299.461182675 -1747.707292580 -1517.034666102
125.710000000
203 40.600000 230.042654659
299.482136309 -1747.682577807 -1517.034195266
125.570000000
###########################################################################################
MY INPUT:
@SET NSTEPS 500000
@SET SOLNAME MOL69
@SET QMCHARGE -1
@SET QMPOISSON MT
@SET FREQ 1
@INCLUDE boxes.inc
&GLOBAL
PRINT_LEVEL LOW
PROJECT system
PROGRAM CP2K
RUN_TYPE MD
WALLTIME 21000
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS ${NSTEPS}
TIMESTEP 0.20
TEMPERATURE 300.0
&THERMOSTAT
TYPE CSVR
REGION DEFINED
&DEFINE_REGION
MM_SUBSYS ATOMIC
&END DEFINE_REGION
&DEFINE_REGION
QM_SUBSYS ATOMIC
&END DEFINE_REGION
&CSVR
TIMECON 50
&END CSVR
&END THERMOSTAT
&PRINT
FORCE_LAST T
&ENERGY
&EACH
MD 1
&END EACH
FILENAME =system.ene
&END ENERGY
&END PRINT
&END MD
&PRINT
&TRAJECTORY SILENT
&EACH
MD ${FREQ}
&END EACH
FILENAME =system.dcd
FORMAT DCD
&END TRAJECTORY
&VELOCITIES OFF
&EACH
MD ${FREQ}
&END EACH
FILENAME =system.vel
&END VELOCITIES
&FORCES
&EACH
MD ${FREQ}
&END EACH
FILENAME =system.for.DCD
FORMAT DCD
&END FORCES
&RESTART_HISTORY OFF
&END RESTART_HISTORY
&RESTART SILENT
ADD_LAST NUMERIC
&EACH
MD ${FREQ}
&END EACH
FILENAME =system.restart
&END RESTART
&END PRINT
&CONSTRAINT
&G3X3
MOLNAME ${SOLNAME}
DISTANCES 1.7958 1.7958 2.85444
ATOMS 1 2 3
EXCLUDE_QM
&END G3X3
# @INCLUDE ./fix_qm_heavy.inc
&END CONSTRAINT
# @INCLUDE ./FREE_ENERGY.inc
&END MOTION
&FORCE_EVAL
METHOD QMMM
&DFT
BASIS_SET_FILE_NAME ./BASIS_SET
POTENTIAL_FILE_NAME ./GTH_POTENTIALS
CHARGE ${QMCHARGE}
&SCF
MAX_SCF 40
EPS_SCF 1.0E-7
SCF_GUESS ATOMIC
&OUTER_SCF
EPS_SCF 1.0E-7
MAX_SCF 10
&END OUTER_SCF
&OT
PRECONDITIONER FULL_ALL
MINIMIZER DIIS
N_DIIS 7
&END OT
&PRINT
&RESTART
ADD_LAST NUMERIC
&EACH
MD ${FREQ}
&END EACH
FILENAME =RESTART
&END RESTART
&RESTART_HISTORY OFF
&END RESTART_HISTORY
&END PRINT
&END SCF
&QS
EPS_DEFAULT 1.0E-12
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
&END QS
&MGRID
CUTOFF 320
COMMENSURATE
&END MGRID
&POISSON
POISSON_SOLVER ${QMPOISSON}
PERIODIC NONE
&END POISSON
&XC
&XC_GRID
XC_SMOOTH_RHO NN50
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&MM
&FORCEFIELD
PARMTYPE AMBER
PARM_FILE_NAME ./molbox.prmtop
EI_SCALE14 0.8333
VDW_SCALE14 0.5
&SPLINE
RCUT_NB 12.0
&END SPLINE
&BOND
KIND AMBER
ATOMS ZN SH
K [angstrom^-2kcalmol] 68.05
R0 [angstrom] 2.311
&END BOND
&BOND
KIND AMBER
ATOMS ZN NB
K [angstrom^-2kcalmol] 45.58
R0 [angstrom] 1.957
&END BOND
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE SPME
ALPHA .35
GMAX ${EWALDA} ${EWALDB} ${EWALDC}
&END EWALD
&END POISSON
&END MM
&QMMM
USE_GEEP_LIB 7
ECOUPL GAUSS
&CELL
ABC ${QMBOXA} ${QMBOXB} ${QMBOXC}
PERIODIC NONE
&END CELL
@INCLUDE ./QM_KIND-LINK.inc
&WALLS
WALL_SKIN [angstrom] 2.0
TYPE REFLECTIVE
&END WALLS
&PRINT
&PROGRAM_RUN_INFO SILENT
&END PROGRAM_RUN_INFO
&PERIODIC_INFO SILENT
&END PERIODIC_INFO
&QMMM_LINK_INFO SILENT
&END QMMM_LINK_INFO
&END PRINT
&END QMMM
&SUBSYS
&CELL
ABC ${MMBOXA} ${MMBOXB} ${MMBOXC}
PERIODIC XYZ
&END CELL
&TOPOLOGY
&MOL_SET
&MERGE_MOLECULES
&BONDS
1520 8881
11567 16020
&END BONDS
&END MERGE_MOLECULES
&END MOL_SET
&GENERATE
&BOND ADD
ATOMS 1520 8881
&END BOND
&BOND ADD
ATOMS 1476 8881
&END BOND
&BOND ADD
ATOMS 961 8881
&END BOND
&BOND ADD
ATOMS 905 8881
&END BOND
&BOND ADD
ATOMS 11567 16020
&END BOND
&BOND ADD
ATOMS 11611 16020
&END BOND
&BOND ADD
ATOMS 11052 16020
&END BOND
&BOND ADD
ATOMS 10996 16020
&END BOND
&END GENERATE
# COORD_FILE_NAME ./molbox.rst7
# COORD_FILE_FORMAT CRD
CONN_FILE_FORMAT psf
CONN_FILE_NAME ./molbox_mod.psf
MOL_CHECK T
&END TOPOLOGY
# @INCLUDE ./COLVARS.inc
@INCLUDE ./KIND_BASIS_POTENTIAL.inc
&END SUBSYS
&END FORCE_EVAL
&EXT_RESTART
RESTART_FILE_NAME system.restart
RESTART_DEFAULT F
RESTART_POS T
RESTART_QMMM F
RESTART_CELL F
RESTART_VEL T
RESTART_CONSTRAINT F
&END EXT_RESTART
one step
MD_ENERGIES| Initialization proceeding
******************************** GO CP2K GO!
**********************************
INITIAL POTENTIAL ENERGY[hartree] =
-0.174799785808E+04
TOTAL INITIAL KINETIC ENERGY[hartree] =
0.230415023002E+03
QM INITIAL KINETIC ENERGY[hartree] =
0.158480716859E+00
TOTAL INITIAL TEMPERATURE[K]
= 299.967
QM INITIAL TEMPERATURE[K]
= 383.480
INITIAL VOLUME[bohr^3] =
0.162855828148E+08
INITIAL CELL LNTHS[bohr] = 0.2427467E+03 0.1883698E+03
0.3561548E+03
INITIAL CELL ANGLS[deg] = 0.9000000E+02 0.9000000E+02
0.9000000E+02
******************************** GO CP2K GO!
**********************************
Translating the system in order to center the QM fragment in the QM
box.
QMMM| Information on the QM/MM Electrostatic Potential:
QMMM| QM/MM Coupling computed collocating the Gaussian Potential
Functions.
ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.):
-1107.691789517704819
Number of
electrons: 282
Number of occupied
orbitals: 141
Number of molecular
orbitals: 141
Number of orbital
functions: 853
Number of independent orbital
functions: 853
Parameters for the always stable predictor-corrector (ASPC) method:
ASPC order: 0
B(1) = 1.000000
Extrapolation method: ASPC
SCF WAVEFUNCTION OPTIMIZATION
----------------------------------- OT
---------------------------------------
----------------------------------- OT
---------------------------------------
Allowing for rotations: F
Optimizing orbital energies: F
Minimizer : DIIS : direct inversion
in the iterative subspace
using : - 7 DIIS vectors
- safer DIIS on
Preconditioner : FULL_ALL : diagonalization, state
selective
Precond_solver : DEFAULT
stepsize : 0.15000000
energy_gap : 0.20000000
eps_taylor : 0.10000E-15
max_taylor : 4
mixed_precision : F
----------------------------------- OT
---------------------------------------
Step Update method Time Convergence Total
energy Change
------------------------------------------------------------------------------
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
1 OT DIIS 0.15E+00 3.6 0.00013824 -660.0889519260
-6.60E+02
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
2 OT DIIS 0.15E+00 3.5 0.00006528 -660.0891544242
-2.02E-04
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
3 OT DIIS 0.15E+00 3.5 0.00004102 -660.0891852777
-3.09E-05
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
4 OT DIIS 0.15E+00 3.5 0.00000905 -660.0891968633
-1.16E-05
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
5 OT DIIS 0.15E+00 3.5 0.00000363 -660.0891972223
-3.59E-07
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
6 OT DIIS 0.15E+00 3.5 0.00000128 -660.0891972771
-5.49E-08
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
7 OT DIIS 0.15E+00 3.5 0.00000050 -660.0891972836
-6.48E-09
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
8 OT DIIS 0.15E+00 3.5 0.00000017 -660.0891972847
-1.03E-09
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
9 OT DIIS 0.15E+00 3.5 0.00000006 -660.0891972848
-1.36E-10
*** SCF run converged in 9 steps ***
Electronic density on regular grids: -282.0000853101
-0.0000853101
Core density on regular grids: 280.9999999999
-0.0000000001
Total charge density on r-space grids: -1.0000853102
Total charge density g-space grids: -1.0000853102
Overlap energy of the core charge distribution:
0.00000627667060
Self energy of the core charge distribution:
-1589.05360281114645
Core Hamiltonian energy:
468.03211746420584
Hartree energy:
628.55173480393682
Exchange-correlation energy:
-147.56475324920322
QM/MM Electrostatic energy:
-20.0546997693
Total energy:
-660.08919728479384
outer SCF iter = 1 RMS gradient = 0.60E-07 energy =
-660.0891972848
outer SCF loop converged in 1 iterations or 9 steps
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):
-660.089197284806573
QMMM| Evaluating forces on MM atoms due to the:
QMMM| - QM/MM Coupling computed collocating the Gaussian Potential
Functions.
QMMM| QM/MM Nuclear Electrostatic Potential :
19.765411770
QMMM| QMMM Total Energy (QM + QMMM electronic + QMMM nuclear):
-640.323785514
QMMM| MM energy NOT included in the above term! Check for:
FORCE_EVAL ( QMMM )
QMMM| that includes both QM, QMMM and MM energy terms!
ENERGY| Total FORCE_EVAL ( QMMM ) energy (a.u.):
-1748.015575032073230
*******************************************************************************
ENSEMBLE TYPE
= NVT
STEP NUMBER
= 154
TIME [FS] =
30.800000
CONSERVED QNTY =
-0.151703373148E+04
INSTANTANEOUS AVERAGES
CPU [S] =
232.00 136.50
{E-E0}/{k_b*N_at} = -0.206102875908E+04
-0.204766829357E+04
POTENTIAL ENERGY[hartree] = -0.174801557503E+04
-0.174803150569E+04
TOTAL KINETIC ENERGY[hartree]= 0.230424151387E+03
0.230737291553E+03
QM KINETIC ENERGY[hartree] = 0.161605341266E+00
0.131303426565E+00
TOTAL TEMPERATURE[K] =
299.979 300.386
QM TEMPERATURE[K] =
391.041 317.719
*******************************************************************************
Translating the system in order to center the QM fragment in the QM
box.
QMMM| Information on the QM/MM Electrostatic Potential:
QMMM| QM/MM Coupling computed collocating the Gaussian Potential
Functions.
ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.):
-1107.709038972358030
Number of
electrons: 282
Number of occupied
orbitals: 141
Number of molecular
orbitals: 141
Number of orbital
functions: 853
Number of independent orbital
functions: 853
Parameters for the always stable predictor-corrector (ASPC) method:
ASPC order: 0
B(1) = 2.000000
B(2) = -1.000000
Extrapolation method: ASPC
SCF WAVEFUNCTION OPTIMIZATION
----------------------------------- OT
---------------------------------------
Allowing for rotations: F
Optimizing orbital energies: F
Minimizer : DIIS : direct inversion
in the iterative subspace
using : - 7 DIIS vectors
- safer DIIS on
Preconditioner : FULL_ALL : diagonalization, state
selective
Precond_solver : DEFAULT
stepsize : 0.15000000
energy_gap : 0.20000000
eps_taylor : 0.10000E-15
max_taylor : 4
mixed_precision : F
----------------------------------- OT
---------------------------------------
Step Update method Time Convergence Total
energy Change
------------------------------------------------------------------------------
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
1 OT DIIS 0.15E+00 3.6 0.00002521 -660.1535199210
-6.60E+02
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
2 OT DIIS 0.15E+00 3.5 0.00001373 -660.1535252226
-5.30E-06
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
3 OT DIIS 0.15E+00 3.5 0.00001003 -660.1535257786
-5.56E-07
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
4 OT DIIS 0.15E+00 3.5 0.00000223 -660.1535261910
-4.12E-07
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
5 OT DIIS 0.15E+00 3.5 0.00000109 -660.1535262125
-2.14E-08
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
6 OT DIIS 0.15E+00 3.5 0.00000042 -660.1535262171
-4.65E-09
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
7 OT DIIS 0.15E+00 3.5 0.00000013 -660.1535262175
-3.66E-10
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
8 OT DIIS 0.15E+00 3.5 0.00000005 -660.1535262175
-1.77E-11
*** SCF run converged in 8 steps ***
Electronic density on regular grids: -282.0001526135
-0.0001526135
Core density on regular grids: 280.9999999999
-0.0000000001
Total charge density on r-space grids: -1.0001526136
Total charge density g-space grids: -1.0001526136
Overlap energy of the core charge distribution:
0.00000641065336
Self energy of the core charge distribution:
-1589.05360281114645
Core Hamiltonian energy:
468.06681187909601
Hartree energy:
628.52418689557771
Exchange-correlation energy:
-147.57280797171762
QM/MM Electrostatic energy:
-20.1181206200
Total energy:
-660.15352621750253
outer SCF iter = 1 RMS gradient = 0.50E-07 energy =
-660.1535262175
outer SCF loop converged in 1 iterations or 8 steps
Adding QM/MM electrostatic potential to the Kohn-Sham potential.
ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):
-660.153526217483318
QMMM| Evaluating forces on MM atoms due to the:
QMMM| - QM/MM Coupling computed collocating the Gaussian Potential
Functions.
QMMM| QM/MM Nuclear Electrostatic Potential :
19.828156874
QMMM| QMMM Total Energy (QM + QMMM electronic + QMMM nuclear):
-640.325369343
QMMM| MM energy NOT included in the above term! Check for:
FORCE_EVAL ( QMMM )
QMMM| that includes both QM, QMMM and MM energy terms!
ENERGY| Total FORCE_EVAL ( QMMM ) energy (a.u.):
-1748.034408315605106
*******************************************************************************
ENSEMBLE TYPE
= NVT
STEP NUMBER
= 155
TIME [FS] =
31.000000
CONSERVED QNTY =
-0.151703386645E+04
INSTANTANEOUS AVERAGES
CPU [S] =
125.98 136.43
{E-E0}/{k_b*N_at} = -0.206102894244E+04
-0.204775449130E+04
POTENTIAL ENERGY[hartree] = -0.174803440832E+04
-0.174803152442E+04
TOTAL KINETIC ENERGY[hartree]= 0.230383030154E+03
0.230735005995E+03
QM KINETIC ENERGY[hartree] = 0.162486103679E+00
0.131504605127E+00
TOTAL TEMPERATURE[K] =
299.925 300.383
QM TEMPERATURE[K] =
393.172 318.206
*******************************************************************************
Translating the system in order to center the QM fragment in the QM
box.
QMMM| Information on the QM/MM Electrostatic Potential:
QMMM| QM/MM Coupling computed collocating the Gaussian Potential
Functions.
ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.):
-1107.728701366057976
More information about the CP2K-user
mailing list