Cons Qty[a.u.] drift? kinetic + potential = Cons Qty ?
marco.stenta
marco.... at gmail.com
Thu Aug 4 10:24:06 UTC 2011
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?
&XC_GRID
XC_SMOOTH_RHO NN50
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
Do you suggest removing it altogether?
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
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
More information about the CP2K-user
mailing list