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