[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