[CP2K:3398] Cons Qty[a.u.] drift? kinetic + potential = Cons Qty ?
Teodoro Laino
teodor... at gmail.com
Thu Aug 4 07:29:58 UTC 2011
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.
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
On Aug 4, 2011, at 9:06 AM, marco.stenta wrote:
> 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
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> --
> You received this message because you are subscribed to the Google Groups "cp2k" group.
> To post to this group, send email to cp... at googlegroups.com.
> To unsubscribe from this group, send email to cp2k+uns... at googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/cp2k?hl=en.
>
More information about the CP2K-user
mailing list