[CP2K:3400] Re: Cons Qty[a.u.] drift? kinetic + potential = Cons Qty ?

Teodoro Laino teodor... at gmail.com
Thu Aug 4 10:43:23 UTC 2011


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/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

> Any suggestion based on your experience is welcome
> 
> thank you very much in advance
> marco
> 
> 
> 
> On 4 Ago, 09:29, Teodoro Laino <teodor... at gmail.com> wrote:
>> 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                
>> 
>> ...
>> 
>> leggi tutto
> 
> -- 
> 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