[CP2K:3402] Re: Cons Qty[a.u.] drift? kinetic + potential = Cons Qty ?
Teodoro Laino
teodor... at gmail.com
Thu Aug 4 14:56:09 UTC 2011
I think that it is more important what people have written about this topic, rather than spending time in expressing personal opinion.
Lots of nice readings out there.
I can recommend this paper:
http://pubs.acs.org/doi/pdf/10.1021/ct1005455
It's about "Long-Range Electrostatic Effects in QM/MM Studies of Enzymatic Reactions: Application of the Solvated Macromolecule Boundary Potential "
This is a very good work by Thiel and coworkers of 2011.
Another good one is:
http://pubs.acs.org/doi/abs/10.1021/ct049941i
It's about "An Efficient Linear-Scaling Ewald Method for Long-Range Electrostatic Interactions in Combined QM/MM Calculations"
by Darren (York) and coworkers of 2006.
Starting from these papers you can find extremely useful information in the therein cited references.
Teo
On Aug 4, 2011, at 2:03 PM, marco.stenta wrote:
> OK:
> the periodic stuff was another question, unrelated to the box
> recenter.
> My question was whether a (non periodic QM)/MM simulation is
> acceptable for publication nowadays.
> In that case I can put some effort and try to learn how to setup a
> full periodic system, if it is really worth.
> I want to avoid being unable to reply to a referee criticizing my
> 1milion hour - long simulation because it is not fully periodic.
> I've seen in the past non periodic was common, but what now? CP2K seem
> a pretty powerful tool and I would know your personal educated opinion
> about that: is a (non periodic QM)/MM reaction study (TI Metadynamics)
> prone to criticism for that or not?
>
> thanks a lot
>
> m
>
>
>
> On 4 Ago, 12:43, Teodoro Laino <teodor... at gmail.com> wrote:
>> 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/dfdde2aa4b6...
>>> 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
>>
>> ...
>>
>> 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.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20110804/48d362f2/attachment.htm>
More information about the CP2K-user
mailing list