Cons Qty[a.u.] drift? kinetic + potential = Cons Qty ?
2,035 views
Skip to first unread message
marco.stenta
Aug 4, 2011, 3:06:38 AM8/4/11
to cp2k
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
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
FORCE_LAST T
&ENERGY
&EACH
MD 1
&END EACH
FILENAME =system.ene
&END ENERGY
&END PRINT
&END MD
&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
&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
&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
Teodoro Laino
Aug 4, 2011, 3:29:58 AM8/4/11
to cp...@googlegroups.com
My hypothesis: The very basic answer is that your QM system is continuously translated in the middle of the QM box.
You say that you keep all the QM atom frozen. In the input file you've sent I don't see that (i.e. it is commented, modulo my personal mistakes in looking at an input of 300 lines!). So I do not think that you are really keeping all QM atoms frozen. This translation is definitely accounting for a constant drift (it's quite normal to observe in PW based codes).
You can try to disable the translation and see what happens.
You say that you keep all the QM atom frozen. In the input file you've sent I don't see that (i.e. it is commented, modulo my personal mistakes in looking at an input of 300 lines!). So I do not think that you are really keeping all QM atoms frozen. This translation is definitely accounting for a constant drift (it's quite normal to observe in PW based codes).
You can try to disable the translation and see what happens.
XC_SMOOTH_RHO may also lead to a non-perfect energy conservation. Why are you using that?
There could be of course other reasons.. but.. let's proceed stepwise..
Regarding the cons. Qty, It should be clear that in an NVT ensemble (compared to NVE) the conserved quantity is not only the sum of the kinetic and potential energy of the particles, but you need to include as well the thermostat energies (which are not printed in the ener file)
The Cons Qty (i.e. including the thermostat energies) is the quantity which should be conserved.
As a side comment: if you're curious about your energy conservation just run a simple NVE. In this case you are summing too many effects that (I kind of reckon) you can't control very well.
Regards,
Teo
> --
> You received this message because you are subscribed to the Google Groups "cp2k" group.
> To post to this group, send email to cp...@googlegroups.com.
> To unsubscribe from this group, send email to cp2k+uns...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/cp2k?hl=en.
>
marco.stenta
Aug 4, 2011, 6:24:06 AM8/4/11
to cp2k
Dear Teodoro,
thanks for the quick answer.
in this simulation the heavy QM atoms are free> i had them frozen in
the previous 2 ps of QMMM-MD to enforce a smooth transition between
the MM-MD and the QMMM-MD: now all is free.
I am using the MT poisson sover, so I understand it is safe to avoid
recentering the QM box as it is translational invariant (CP2K_INPUT /
FORCE_EVAL / QMMM / NOCENTER NOCENTER0, if I am right): I'll do that
I am using XC_SMOOTH_RHO because I've seen it used often (It is just
guess work for some keyword) often XC_SMOOTH_RHO NN10,
but I've also seen this post by Hutter where he suggest using it
https://groups.google.com/group/cp2k/browse_thread/thread/75ce24ccda8dc3ef/4444267819ea4acd?hl=it&lnk=gst&q=conserved+quantity#4444267819ea4acd
quote(
related to the numerics of the BLYP functional. You can either
increase the energy cutoff from 280 Ry to a considerable higher
value or use a smoothing procedure of the density.
XC_SMOOTH_RHO NN50
would be my choice.
)unquote
but I've seen you suggest instead leaving the defaults (XC_DERIV PW)
https://groups.google.com/group/cp2k/browse_thread/thread/bb9d011f11d5cf31/d4fdff84d4394fa1?hl=it&lnk=gst&q=XC_SMOOTH_RHO#d4fdff84d4394fa1
so what should I do with the XC_grid section?
what about EPSFIT?
now the defaults is 1.00000000E-04
but here it was suggested to use 10^-2
https://groups.google.com/group/cp2k/browse_thread/thread/dfdde2aa4b64c8f/40181675b6896106?hl=it&lnk=gst&q=conserved+quantity#40181675b6896106
which is the most accurate?
Then another point:
I'm familiar with QMMM where the QM is not periodic.
Do you suggest me to switch to a periodic description? I is a standard
enzymatic system on which I will study a reaction mechanism with TI
and metadynamics: are the results obtained with a (non-periodic QM)/
MM simulation acceptable?
What do you usually do in such cases?
I'll definitely check with an NVE for the energy conservation.
but how can I assess later equilibration in my NVT system if the
conserved quantity still drifts? is the sum of kinetic and potential a
good parameter? RMSD of the forces of the active site?
Is conserved quantity a parameter to worry about?
Any suggestion based on your experience is welcome
thank you very much in advance
marco
> ...
>
> leggi tutto
thanks for the quick answer.
in this simulation the heavy QM atoms are free> i had them frozen in
the previous 2 ps of QMMM-MD to enforce a smooth transition between
the MM-MD and the QMMM-MD: now all is free.
I am using the MT poisson sover, so I understand it is safe to avoid
recentering the QM box as it is translational invariant (CP2K_INPUT /
FORCE_EVAL / QMMM / NOCENTER NOCENTER0, if I am right): I'll do that
I am using XC_SMOOTH_RHO because I've seen it used often (It is just
guess work for some keyword) often XC_SMOOTH_RHO NN10,
but I've also seen this post by Hutter where he suggest using it
https://groups.google.com/group/cp2k/browse_thread/thread/75ce24ccda8dc3ef/4444267819ea4acd?hl=it&lnk=gst&q=conserved+quantity#4444267819ea4acd
quote(
related to the numerics of the BLYP functional. You can either
increase the energy cutoff from 280 Ry to a considerable higher
value or use a smoothing procedure of the density.
XC_SMOOTH_RHO NN50
would be my choice.
)unquote
but I've seen you suggest instead leaving the defaults (XC_DERIV PW)
https://groups.google.com/group/cp2k/browse_thread/thread/bb9d011f11d5cf31/d4fdff84d4394fa1?hl=it&lnk=gst&q=XC_SMOOTH_RHO#d4fdff84d4394fa1
so what should I do with the XC_grid section?
&XC_GRID
XC_SMOOTH_RHO NN50
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
Do you suggest removing it altogether?
XC_SMOOTH_RHO NN50
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
what about EPSFIT?
now the defaults is 1.00000000E-04
but here it was suggested to use 10^-2
https://groups.google.com/group/cp2k/browse_thread/thread/dfdde2aa4b64c8f/40181675b6896106?hl=it&lnk=gst&q=conserved+quantity#40181675b6896106
which is the most accurate?
Then another point:
I'm familiar with QMMM where the QM is not periodic.
Do you suggest me to switch to a periodic description? I is a standard
enzymatic system on which I will study a reaction mechanism with TI
and metadynamics: are the results obtained with a (non-periodic QM)/
MM simulation acceptable?
What do you usually do in such cases?
I'll definitely check with an NVE for the energy conservation.
but how can I assess later equilibration in my NVT system if the
conserved quantity still drifts? is the sum of kinetic and potential a
good parameter? RMSD of the forces of the active site?
Is conserved quantity a parameter to worry about?
Any suggestion based on your experience is welcome
thank you very much in advance
marco
>
> leggi tutto
Teodoro Laino
Aug 4, 2011, 6:43:23 AM8/4/11
to cp...@googlegroups.com
We provide suggestion: finally you need to do things at your own convenience.
Having said that:
Having said that:
>
> so what should I do with the XC_grid section?
> &XC_GRID
> XC_SMOOTH_RHO NN50
> XC_DERIV SPLINE2_SMOOTH
> &END XC_GRID
>
> Do you suggest removing it altogether ?
It's your decision: you found two different posts, analyze them, do tests (!) and draw your conclusions.
>
> what about EPSFIT?
> now the defaults is 1.00000000E-04
> but here it was suggested to use 10^-2
>
> https://groups.google.com/group/cp2k/browse_thread/thread/dfdde2aa4b64c8f/40181675b6896106?hl=it&lnk=gst&q=conserved+quantity#40181675b6896106
> which is the most accurate?
EPSFIT regards GAPW only. You're using GPW. Why are interested in that?
> Then another point:
> I'm familiar with QMMM where the QM is not periodic.
> Do you suggest me to switch to a periodic description? I is a standard
> enzymatic system on which I will study a reaction mechanism with TI
> and metadynamics: are the results obtained with a (non-periodic QM)/
> MM simulation acceptable?
> What do you usually do in such cases?
I've never talked in my reply (to your message) about periodic or non-periodic. I've only said that you should try to:
-) switch off the translation/centering of the QM system
-) use an NVE to inspect the energy conservation
the translation/centering of the QM system does not depend on the periodic or non-periodic. It depends on the fact that the density is contained in a simulation box and cannot travel freely around space.
Moreover, since you have QM atoms moving, the center of QM system will be constantly moved in the center of the QM box and this creates the drift.
Periodic / non-periodic is a totally different story and is not affecting the conserved qty with a constant drift (maybe in another way).
> I'll definitely check with an NVE for the energy conservation.
> but how can I assess later equilibration in my NVT system if the
> conserved quantity still drifts? is the sum of kinetic and potential a
> good parameter? RMSD of the forces of the active site?
> Is conserved quantity a parameter to worry about?
What about the conserved quantity printed in the .ener file (second-last column)? That contains everything..
If it drifts there is a problem with your setup. In this case, as I said, very probably is the continuous translation/centering of the QM system.
Teo
marco.stenta
Aug 4, 2011, 8:03:35 AM8/4/11
to cp2k
OK:
the periodic stuff was another question, unrelated to the box
recenter.
My question was whether a (non periodic QM)/MM simulation is
acceptable for publication nowadays.
In that case I can put some effort and try to learn how to setup a
full periodic system, if it is really worth.
I want to avoid being unable to reply to a referee criticizing my
1milion hour - long simulation because it is not fully periodic.
I've seen in the past non periodic was common, but what now? CP2K seem
a pretty powerful tool and I would know your personal educated opinion
about that: is a (non periodic QM)/MM reaction study (TI Metadynamics)
prone to criticism for that or not?
thanks a lot
m
On 4 Ago, 12:43, Teodoro Laino <teodoro.la...@gmail.com> wrote:
> We provide suggestion: finally you need to do things at your own convenience.
> Having said that:
>
>
>
> > so what should I do with the XC_grid section?
> > &XC_GRID
> > XC_SMOOTH_RHO NN50
> > XC_DERIV SPLINE2_SMOOTH
> > &END XC_GRID
>
> > Do you suggest removing it altogether ?
>
> It's your decision: you found two different posts, analyze them, do tests (!) and draw your conclusions.
>
>
>
> > what about EPSFIT?
> > now the defaults is 1.00000000E-04
> > but here it was suggested to use 10^-2
>
> >https://groups.google.com/group/cp2k/browse_thread/thread/dfdde2aa4b6...
> ...
>
> leggi tutto
the periodic stuff was another question, unrelated to the box
recenter.
My question was whether a (non periodic QM)/MM simulation is
acceptable for publication nowadays.
In that case I can put some effort and try to learn how to setup a
full periodic system, if it is really worth.
I want to avoid being unable to reply to a referee criticizing my
1milion hour - long simulation because it is not fully periodic.
I've seen in the past non periodic was common, but what now? CP2K seem
a pretty powerful tool and I would know your personal educated opinion
about that: is a (non periodic QM)/MM reaction study (TI Metadynamics)
prone to criticism for that or not?
thanks a lot
m
On 4 Ago, 12:43, Teodoro Laino <teodoro.la...@gmail.com> wrote:
> We provide suggestion: finally you need to do things at your own convenience.
> Having said that:
>
>
>
> > so what should I do with the XC_grid section?
> > &XC_GRID
> > XC_SMOOTH_RHO NN50
> > XC_DERIV SPLINE2_SMOOTH
> > &END XC_GRID
>
> > Do you suggest removing it altogether ?
>
> It's your decision: you found two different posts, analyze them, do tests (!) and draw your conclusions.
>
>
>
> > what about EPSFIT?
> > now the defaults is 1.00000000E-04
> > but here it was suggested to use 10^-2
>
>
> leggi tutto
Teodoro Laino
Aug 4, 2011, 10:56:09 AM8/4/11
to cp...@googlegroups.com
I think that it is more important what people have written about this topic, rather than spending time in expressing personal opinion.
Lots of nice readings out there.
I can recommend this paper:
It's about "Long-Range Electrostatic Effects in QM/MM Studies of Enzymatic Reactions: Application of the Solvated Macromolecule Boundary Potential "
This is a very good work by Thiel and coworkers of 2011.
Another good one is:
It's about "An Efficient Linear-Scaling Ewald Method for Long-Range Electrostatic Interactions in Combined QM/MM Calculations"
by Darren (York) and coworkers of 2006.
Starting from these papers you can find extremely useful information in the therein cited references.
Teo
marco.stenta
Aug 4, 2011, 11:05:16 AM8/4/11
to cp2k
Thanks,
I'll give a look to them,
cheers
m
> ...
>
> leggi tutto
I'll give a look to them,
cheers
m
>
> leggi tutto
Reply all
Reply to author
Forward
0 new messages