Cons. Qty for different ensembles + Pressure fluctuation

[Léon Luntadila Lufungula]

Dear all,

I ran several simulations on a box of 32 water molecules with the TIP3P force field:

1) 1 ns NVE trajectory with timestep 0.5 fs

2) 1 ns NVT trajectory at 298.15K with timestep 0.5 fs and thermostat time constant of 100fs

 3) 1 ns NPT_I trajectory restarted from the final structure of the NVT trajectory at 298.15K and 1 bar with timestep 0.5 and thermostat and barostat time constants of 100 fs

I plotted the Cons. Qty [a.u.] vs. Time [fs] for the three simulations and saw that the constant of motion was only more or less constant for the NVE trajectory. The NVT and NPT trajectories both had strong oscillations and significant drift of the constant of motion (see below) and I was wondering if this is normal behaviour or that this means that the system is not yet equilibrated? I would think that the Cons. Qty. should remain constant but I can also understand that the thermostat and barostat might influence its behaviour in a way.

NVE:

NVT:

NPT:

The simulations of the water box were small test cases, but I also need to know what behaviour to expect when going to my production runs with my systems of interest...

Additionally, in the NPT ensemble, the pressure seems to fluctuate very strongly from -20000 to +20000 bar which seems to be quite excessive to me. Is this normal? I've also seen this asked in another post, but that person did not get a reply unfortunately....

Thanks in advance for the help!

Kind regards,

Léon

[Krack Matthias]

Dear Léon

 

a 32 H2O cell is still a relatively small system which can cause large fluctuations. I suggest to try at least the 6x6x6 = 216 H2O input. The sample should be equilibrated well in advance before starting the sampling run.

 

HTH

 

Matthias

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Monday, 22 January 2024 at 12:56
To: cp2k <cp...@googlegroups.com>
Subject: [CP2K:19784] Cons. Qty for different ensembles + Pressure fluctuation

Dear all,

 

I ran several simulations on a box of 32 water molecules with the TIP3P force field:

1) 1 ns NVE trajectory with timestep 0.5 fs

2) 1 ns NVT trajectory at 298.15K with timestep 0.5 fs and thermostat time constant of 100fs

 3) 1 ns NPT_I trajectory restarted from the final structure of the NVT trajectory at 298.15K and 1 bar with timestep 0.5 and thermostat and barostat time constants of 100 fs

 

I plotted the Cons. Qty [a.u.] vs. Time [fs] for the three simulations and saw that the constant of motion was only more or less constant for the NVE trajectory. The NVT and NPT trajectories both had strong oscillations and significant drift of the constant of motion (see below) and I was wondering if this is normal behaviour or that this means that the system is not yet equilibrated? I would think that the Cons. Qty. should remain constant but I can also understand that the thermostat and barostat might influence its behaviour in a way.

NVE:

NVT:

NPT:

The simulations of the water box were small test cases, but I also need to know what behaviour to expect when going to my production runs with my systems of interest...

 

Additionally, in the NPT ensemble, the pressure seems to fluctuate very strongly from -20000 to +20000 bar which seems to be quite excessive to me. Is this normal? I've also seen this asked in another post, but that person did not get a reply unfortunately....

 

Thanks in advance for the help!

 

Kind regards,

Léon

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/76d3e94b-4a67-4fcc-8861-78a566ebb107n%40googlegroups.com.

[Léon Luntadila Lufungula]

Dear Matthias,

I have similar results for a box of 256 water molecules:

NVT:

NPT: (fluctuations are a bit smaller but still in the range of several kbar...)

You say that the sample should be equilibrated well in advance before starting the sampling run, so do you mean I should run an NVT simulation for longer until it equilibrates before running an NPT from the thermalized end structure? Or do you mean both NVT and NPT should be run for longer until the Cons. Qty. is also equilibrated? Because in both cases the temperature seems to be equilibrated quite well and the pressure in the NPT ensemble is also constant albeit with very large oscillations...

Kind regards,

Léon

[Léon Luntadila Lufungula]

P.S. These calculations were done based of the input files in https://www.cp2k.org/exercises:2018_uzh_acpc2:mol_sol, because the Fist regtest-water folder was not yet available at the time, but the FF params are correct, the coordinates are just different than the one molecule box in the tests directory.

[Jürg Hutter]

Hi

1) The fluctuations of the pressure have been reduced. This is expected
as the fluctuations should be proportional to the sqrt of system size.

2) What thermostats/barostats are you using and what are the corresponding
settings?

regards
JH

________________________________________

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>

Sent: Monday, January 22, 2024 2:03 PM
To: cp2k
Subject: Re: [CP2K:19787] Cons. Qty for different ensembles + Pressure fluctuation

Dear Matthias,

I have similar results for a box of 256 water molecules:
NVT:

[NVT-256.png]

NPT: (fluctuations are a bit smaller but still in the range of several kbar...)

[NPT-256.png][Pressure-256.png]

You say that the sample should be equilibrated well in advance before starting the sampling run, so do you mean I should run an NVT simulation for longer until it equilibrates before running an NPT from the thermalized end structure? Or do you mean both NVT and NPT should be run for longer until the Cons. Qty. is also equilibrated? Because in both cases the temperature seems to be equilibrated quite well and the pressure in the NPT ensemble is also constant albeit with very large oscillations...

Kind regards,
Léon


On Monday 22 January 2024 at 13:40:05 UTC+1 Krack Matthias wrote:
Dear Léon

a 32 H2O cell is still a relatively small system which can cause large fluctuations. I suggest to try at least the 6x6x6 = 216 H2O input<https://github.com/cp2k/cp2k/blob/master/tests/Fist/regtest-water/H2O-TIP3P-6x6x6.inp>. The sample should be equilibrated well in advance before starting the sampling run.

HTH

Matthias

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Monday, 22 January 2024 at 12:56
To: cp2k <cp...@googlegroups.com>
Subject: [CP2K:19784] Cons. Qty for different ensembles + Pressure fluctuation
Dear all,

I ran several simulations on a box of 32 water molecules with the TIP3P force field:
1) 1 ns NVE trajectory with timestep 0.5 fs
2) 1 ns NVT trajectory at 298.15K with timestep 0.5 fs and thermostat time constant of 100fs
3) 1 ns NPT_I trajectory restarted from the final structure of the NVT trajectory at 298.15K and 1 bar with timestep 0.5 and thermostat and barostat time constants of 100 fs

I plotted the Cons. Qty [a.u.] vs. Time [fs] for the three simulations and saw that the constant of motion was only more or less constant for the NVE trajectory. The NVT and NPT trajectories both had strong oscillations and significant drift of the constant of motion (see below) and I was wondering if this is normal behaviour or that this means that the system is not yet equilibrated? I would think that the Cons. Qty. should remain constant but I can also understand that the thermostat and barostat might influence its behaviour in a way.
NVE:

[https://groups.google.com/group/cp2k/attach/720b7b929b503/NVE.png?part=0.1&view=1]
NVT:
[https://groups.google.com/group/cp2k/attach/720b7b929b503/NVT.png?part=0.4&view=1]
NPT:
[https://groups.google.com/group/cp2k/attach/720b7b929b503/NPT.png?part=0.3&view=1]

The simulations of the water box were small test cases, but I also need to know what behaviour to expect when going to my production runs with my systems of interest...

Additionally, in the NPT ensemble, the pressure seems to fluctuate very strongly from -20000 to +20000 bar which seems to be quite excessive to me. Is this normal? I've also seen this asked in another post<https://groups.google.com/g/cp2k/c/gxA8r6wmyF4/m/k8QT7iglAAAJ>, but that person did not get a reply unfortunately....
[https://groups.google.com/group/cp2k/attach/720b7b929b503/Pressure.png?part=0.2&view=1]

Thanks in advance for the help!

Kind regards,
Léon
--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/76d3e94b-4a67-4fcc-8861-78a566ebb107n%40googlegroups.com<https://groups.google.com/d/msgid/cp2k/76d3e94b-4a67-4fcc-8861-78a566ebb107n%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
You received this message because you are subscribed to the Google Groups "cp2k" group.

To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com<mailto:cp2k+uns...@googlegroups.com>.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/26a4ceaf-64b9-4335-aae6-b163bf2aa5ccn%40googlegroups.com<https://groups.google.com/d/msgid/cp2k/26a4ceaf-64b9-4335-aae6-b163bf2aa5ccn%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Krack Matthias]

Hi

 

The GMAX parameter for Ewald is a bit sloppy in these inputs. Suggested values are 1 to 2 points per Angstrom based on the cell edge lengths, i.e. GMAX 20 is numerically more appropriate. In the example I gave you the GMAX value is automatically adapted to the cell size.

 

Best

 

Matthias

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Monday, 22 January 2024 at 14:13
To: cp2k <cp...@googlegroups.com>
Subject: Re: [CP2K:19788] Cons. Qty for different ensembles + Pressure fluctuation

P.S. These calculations were done based of the input files in https://www.cp2k.org/exercises:2018_uzh_acpc2:mol_sol, because the Fist regtest-water folder was not yet available at the time, but the FF params are correct, the coordinates are just different than the one molecule box in the tests directory.

 

On Monday 22 January 2024 at 14:03:06 UTC+1 Léon Luntadila Lufungula wrote:

Dear Matthias,

 

I have similar results for a box of 256 water molecules:

NVT:

NPT: (fluctuations are a bit smaller but still in the range of several kbar...)

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/e6d4002e-91a7-45cd-817c-adc66616d1d6n%40googlegroups.com.

[Léon Luntadila Lufungula]

Hi Jurg and Matthias,

Jurg:

1) Okay, good to know that I can reduce the fluctuations by increasing the system size. Still, are fluctuations in the order of kbar not too large when the desired pressure is 1 bar?

2) I'm using the CSVR thermostat with a timecon of 100 fs and the barostat is the standard one with a timecon of 100 fs by just including the following lines to the input and changing the ensemble to NPT_I: 

&BAROSTAT

  TIMECON 100

&END BAROSTAT

2-continued) I do see now that the BAROSTAT section also has a THERMOSTAT subsection. Should the thermostat be defined here in the case of an NPT calculation or is it in the THERMOSTAT subsection of MD as is the case for an NVT run? If not, what is the use of the BAROSTAT/THERMOSTAT subsection?

Matthias:

Thanks for pointing this out to me! I remembered noticing that the value for GMAX was too small this weekend, but I forgot to adjust it accordingly. I'll try to re-run the calculations with the inputs in the regtest-water directory to see if I get better results! Maybe one more question if that's okay for you. I'm also running / getting ready to run AIMD simulations and these are quite heavy so it takes a long time to simulate even a couple of ps. If I see a drift in the constant of motion like in the pictures above, is there a way of knowing if this is due to an error in the input or just because the system has not equilibrated yet? Or is this just a matter of simulating for longer and longer times to see if it goes towards an equilibrium or keeps drifting? Perhaps the "Energy drift per atom [K]" in the output can be of some help, but I'm not certain how it's calculated or how I can use it to check my trajectory...

Thanks for all the help!

Kind regards,

Léon

[Krack Matthias]

Hi Léon

 

Indeed, the "Energy drift per atom [K]" provides a good indicator for checking the equilibration. Given the small cell sizes used in atomistic simulations, pressure fluctuation up to a kbar range are normal, because a mHa/Bohr^3 corresponds to roughly 294 kbar. The pressure usually equilibrates much more slowly than the temperature. On restart, with RESTART_AVERAGES NO, the averages should evolve smoothly and should not change too much for an equilibrated system.

 

Best

 

Matthias

 

.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/9ee8f63c-639e-4708-8a49-19f2cd8b87e4n%40googlegroups.com.

[Léon Luntadila Lufungula]

Hi Matthias,

1) Okay, I'll keep an eye on the "Energy drift per atom [K]"-values then for my AIMD simulations, I'm guessing for an equilibrated structure these should be near or below 1K? 

2) That's a relief that fluctuations in the order of kbar are normal, I was getting worried about what to do to get them in the order of bar!

3) I'll make sure to set RESTART_AVERAGES to NO when restarting my NPT from my NVT, or is it actually not necessary to first do an NVT before running NPT?

Kind regards,

Léon

[Krack Matthias]

1. Yes, if possible
2. The average pressure wrt the desired target pressure however should be in the order of +/- 20 bar
3. You can start directly with NpT. I would add “RESTART_COUNTERS no” to “RESTART_AVERAGES no”.

 

NPT: (fluctuations are a bit smaller but still in the range of several kbar...)

 

You say that the sample should be equilibrated well in advance before starting the sampling run, so do you mean I should run an NVT simulation for longer until it equilibrates before running an NPT from the thermalized end structure? Or do you mean both NVT and NPT should be run for longer until the Cons. Qty. is also equilibrated? Because in both cases the temperature seems to be equilibrated quite well and the pressure in the NPT ensemble is also constant albeit with very large oscillations...

 

Kind regards,

Léon

 

 

On Monday 22 January 2024 at 13:40:05 UTC+1 Krack Matthias wrote:

Dear Léon

 

a 32 H2O cell is still a relatively small system which can cause large fluctuations. I suggest to try at least the 6x6x6 = 216 H2O input. The sample should be equilibrated well in advance before starting the sampling run.

 

HTH

 

Matthias

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Monday, 22 January 2024 at 12:56
To: cp2k <cp...@googlegroups.com>
Subject: [CP2K:19784] Cons. Qty for different ensembles + Pressure fluctuation

Dear all,

 

I ran several simulations on a box of 32 water molecules with the TIP3P force field:

1) 1 ns NVE trajectory with timestep 0.5 fs

2) 1 ns NVT trajectory at 298.15K with timestep 0.5 fs and thermostat time constant of 100fs

 3) 1 ns NPT_I trajectory restarted from the final structure of the NVT trajectory at 298.15K and 1 bar with timestep 0.5 and thermostat and barostat time constants of 100 fs

 

I plotted the Cons. Qty [a.u.] vs. Time [fs] for the three simulations and saw that the constant of motion was only more or less constant for the NVE trajectory. The NVT and NPT trajectories both had strong oscillations and significant drift of the constant of motion (see below) and I was wondering if this is normal behaviour or that this means that the system is not yet equilibrated? I would think that the Cons. Qty. should remain constant but I can also understand that the thermostat and barostat might influence its behaviour in a way.

NVE:

NVT:

NPT:

The simulations of the water box were small test cases, but I also need to know what behaviour to expect when going to my production runs with my systems of interest...

 

Additionally, in the NPT ensemble, the pressure seems to fluctuate very strongly from -20000 to +20000 bar which seems to be quite excessive to me. Is this normal? I've also seen this asked in another post, but that person did not get a reply unfortunately....

 

Thanks in advance for the help!

 

Kind regards,

Léon

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/76d3e94b-4a67-4fcc-8861-78a566ebb107n%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/e6d4002e-91a7-45cd-817c-adc66616d1d6n%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com

.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/9ee8f63c-639e-4708-8a49-19f2cd8b87e4n%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/e9d20a23-e4e5-4a63-a807-6858c86a9b29n%40googlegroups.com.

[Léon Luntadila Lufungula]

Hi Matthias,

So I started an NpT directly from the inputs you gave me with a box containing 216 water molecules. After a 300ps equilibration (plots thus starrt at 0.3 ns), the pressure fluctuations are again in the 10kbar range which is okay as you pointed out, while the average pressure is around -14 bar, so tha'ts also within the 20 bar range you mentioned. The energy drift per atom is also below 1K so everything seems to work reasonably. 

I'm starting my AIMD simulations as we speak, so hopefully all goes well there too. I read that it is best to do a short equilibration with a strong coupling to the thermostat (e.g. CSVR, TIMECON 20 and REGION MASSIVE), before going to a weaker coupling (e.g. CSVR, TIMECON 200 and REGION GLOBAL), but is there also a rule of thumb w.r.t. the barostat TIMECON? I regularly see TIMECON 1000, but is it useful to set it to TIMECON 100 in this early equilibration?

Thanks for all the help already, it really helps me a lot given the fact that I don't have a senior researcher with a lot of hands-on CP2K experience around!

Kind regards,

Léon

[Krack Matthias]

Dear Léon

 

Using a MASSIVE CSVR thermostat with a small TIMECON of about 50 is appropriate and efficient for equilibration. Thereafter, the TIMECON for the thermostat can be increased to 500-1000 for the sampling (production) runs and also a GLOBAL instead of a MASSIVE thermostat can be used, if possible, to reduce the impact of the thermostat on the dynamics. You will find different opinions about the absolute and relative values for the TIMECONs of thermostat and barostat in the web or literature. I suggest to use a two to five times larger TIMECON for the barostat compared to the thermostat. You have, however, to keep in mind that the appropriate (best) parameters for baro- and thermostat depend on the type of system and the selected conditions eventually.

 

HTH

 

Matthias  

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Tuesday, 23 January 2024 at 11:47
To: cp2k <cp...@googlegroups.com>
Subject: Re: [CP2K:19801] Cons. Qty for different ensembles + Pressure fluctuation

Hi Matthias,

 

So I started an NpT directly from the inputs you gave me with a box containing 216 water molecules. After a 300ps equilibration (plots thus starrt at 0.3 ns), the pressure fluctuations are again in the 10kbar range which is okay as you pointed out, while the average pressure is around -14 bar, so tha'ts also within the 20 bar range you mentioned. The energy drift per atom is also below 1K so everything seems to work reasonably. 

 

I'm starting my AIMD simulations as we speak, so hopefully all goes well there too. I read that it is best to do a short equilibration with a strong coupling to the thermostat (e.g. CSVR, TIMECON 20 and REGION MASSIVE), before going to a weaker coupling (e.g. CSVR, TIMECON 200 and REGION GLOBAL), but is there also a rule of thumb w.r.t. the barostat TIMECON? I regularly see TIMECON 1000, but is it useful to set it to TIMECON 100 in this early equilibration?

[Léon Luntadila Lufungula]

Dear Matthias,

My AIMD calculations are on a slab wherein I keep my bottom layer fixed, but apparently REGION MASSIVE is not compatible with constraints, so I'll have to use REGION GLOBAL. Either way, thanks a lot for the detailed answer! For a novice MD practicioner as myself this is already enough information to get started! I'll try the values you suggested during my production run while keeping an eye on the temperature and pressure during the simulation; I guess when I've chosen my TIMECON too large, this will result in the system drifting away from the equilibrium values.

As always, thanks a lot for all the help!

Kind regards,

Léon

[Léon Luntadila Lufungula]

Dear Matthias,

I've also tried to equilibrate a box of 216 water molecules with DFT and while the temperature and pressure seem to equilibrate quite easily, the constant of motion shows large jumps despite the potential and kinetic energy also being quite constant... If I'm not mistaken, the constant of motion in the NPT ensemble is the energy of the system + the energy of the thermostat + energy of barostat, as shown in this lecture from 2016, so this would suggest that the fluctuations are due to the thermostat or barostat, right? Like the TIP3P NPT run, I used a thermostat timecon of 100 fs and a barostat timecon of 1000 fs and after 21 125 fs (blue vertical line in figures), I restarted the calculation and increased the thermostat timecon to 500 fs. I accidentaly set RESTART_BAROSTAT and RESTART_BAROSTAT_THERMOSTAT to FALSE upon restarting, but I'm unsure if this might have lead to the large jumps in the constant of motion later at t > 28 000 fs... 

Same goes for the NPT_F calculation of my solvated TiO2 slab. I did a 10 ps simulation with thermostat timecon 50 and barostat timecon 100 and everything seems to equilibrate except for the constant of motion which seems to steadily rise. After 10 ps of simulation I restarted the calculation and again accidentaly messed up my input and set RESTART_BAROSTAT and RESTART_BAROSTAT_THERMOSTAT to FALSE but this time I didn't even change the timecons of the thermostat and barostat, but this time I immediately see a significant jump in the constant of motion after which it continues to increase... Any idea how I can fix this so that the constant of motion stays...well...constant?

Kind regards,

Léon

[Léon Luntadila Lufungula]

Dear Matthias,

My apologies, I forgot to add the input files. You can find them in attachment to this message.

Kind regards,

Léon

[Krack Matthias]

Dear Léon

 

the constant of motion mentioned in the lecture should correspond to the conserved quantity (6^th) column in the .ener file. Its value will change when you change simulation parameters including those for the thermostat and/or barostat. While the temperature fluctuations of the system might be equilibrated after 10 ps in an acceptable manner, the cell fluctuations will most likely need much longer equilibration times of at least 30 ps or more. Given the computational costs for AIMD, one usually tends to cut corners which results in a significant deterioration of the sampling quality.

 

Best

 

Matthias

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Monday, 29 January 2024 at 10:38
To: cp2k <cp...@googlegroups.com>
Subject: Re: [CP2K:19838] Cons. Qty for different ensembles + Pressure fluctuation

Dear Matthias,

 

I've also tried to equilibrate a box of 216 water molecules with DFT and while the temperature and pressure seem to equilibrate quite easily, the constant of motion shows large jumps despite the potential and kinetic energy also being quite constant... If I'm not mistaken, the constant of motion in the NPT ensemble is the energy of the system + the energy of the thermostat + energy of barostat, as shown in this lecture from 2016, so this would suggest that the fluctuations are due to the thermostat or barostat, right? Like the TIP3P NPT run, I used a thermostat timecon of 100 fs and a barostat timecon of 1000 fs and after 21 125 fs (blue vertical line in figures), I restarted the calculation and increased the thermostat timecon to 500 fs. I accidentaly set RESTART_BAROSTAT and RESTART_BAROSTAT_THERMOSTAT to FALSE upon restarting, but I'm unsure if this might have lead to the large jumps in the constant of motion later at t > 28 000 fs... 

Same goes for the NPT_F calculation of my solvated TiO2 slab. I did a 10 ps simulation with thermostat timecon 50 and barostat timecon 100 and everything seems to equilibrate except for the constant of motion which seems to steadily rise. After 10 ps of simulation I restarted the calculation and again accidentaly messed up my input and set RESTART_BAROSTAT and RESTART_BAROSTAT_THERMOSTAT to FALSE but this time I didn't even change the timecons of the thermostat and barostat, but this time I immediately see a significant jump in the constant of motion after which it continues to increase... Any idea how I can fix this so that the constant of motion stays...well...constant?

 

Kind regards,

Léon

 

On Tuesday 23 January 2024 at 14:41:46 UTC+1 Léon Luntadila Lufungula wrote:

Dear Matthias,

 

My AIMD calculations are on a slab wherein I keep my bottom layer fixed, but apparently REGION MASSIVE is not compatible with constraints, so I'll have to use REGION GLOBAL. Either way, thanks a lot for the detailed answer! For a novice MD practicioner as myself this is already enough information to get started! I'll try the values you suggested during my production run while keeping an eye on the temperature and pressure during the simulation; I guess when I've chosen my TIMECON too large, this will result in the system drifting away from the equilibrium values.

 

As always, thanks a lot for all the help!

 

Kind regards,

Léon

On Tuesday 23 January 2024 at 14:21:05 UTC+1 Krack Matthias wrote:

Dear Léon

 

Using a MASSIVE CSVR thermostat with a small TIMECON of about 50 is appropriate and efficient for equilibration. Thereafter, the TIMECON for the thermostat can be increased to 500-1000 for the sampling (production) runs and also a GLOBAL instead of a MASSIVE thermostat can be used, if possible, to reduce the impact of the thermostat on the dynamics. You will find different opinions about the absolute and relative values for the TIMECONs of thermostat and barostat in the web or literature. I suggest to use a two to five times larger TIMECON for the barostat compared to the thermostat. You have, however, to keep in mind that the appropriate (best) parameters for baro- and thermostat depend on the type of system and the selected conditions eventually.

 

HTH

 

Matthias  

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Tuesday, 23 January 2024 at 11:47
To: cp2k <cp...@googlegroups.com>
Subject: Re: [CP2K:19801] Cons. Qty for different ensembles + Pressure fluctuation

Hi Matthias,

 

So I started an NpT directly from the inputs you gave me with a box containing 216 water molecules. After a 300ps equilibration (plots thus starrt at 0.3 ns), the pressure fluctuations are again in the 10kbar range which is okay as you pointed out, while the average pressure is around -14 bar, so tha'ts also within the 20 bar range you mentioned. The energy drift per atom is also below 1K so everything seems to work reasonably. 

 

I'm starting my AIMD simulations as we speak, so hopefully all goes well there too. I read that it is best to do a short equilibration with a strong coupling to the thermostat (e.g. CSVR, TIMECON 20 and REGION MASSIVE), before going to a weaker coupling (e.g. CSVR, TIMECON 200 and REGION GLOBAL), but is there also a rule of thumb w.r.t. the barostat TIMECON? I regularly see TIMECON 1000, but is it useful to set it to TIMECON 100 in this early equilibration?

 

Thanks for all the help already, it really helps me a lot given the fact that I don't have a senior researcher with a lot of hands-on CP2K experience around!

 

Kind regards,

Léon

 

--

You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/df0cb86d-d2ad-45e0-a054-1fc7d7e0f789n%40googlegroups.com.

[Krack Matthias]

Dear Léon

 

MD based on DFT compared to classical potentials suffers more from numerical noise in the forces. Therefore, you should reduce EPS_DEFAULT from 1.0E-10 to at least 1.0E-12 in order to reduce that noise. Nevertheless, you will always observe a noticeable energy drift when you run long enough.

 

HTH

 

Matthias

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Monday, 29 January 2024 at 11:03
To: cp2k <cp...@googlegroups.com>
Subject: Re: [CP2K:19840] Cons. Qty for different ensembles + Pressure fluctuation

Dear Matthias,

 

My apologies, I forgot to add the input files. You can find them in attachment to this message.

 

Kind regards,

Léon

On Monday 29 January 2024 at 10:36:30 UTC+1 Léon Luntadila Lufungula wrote:

Dear Matthias,

 

I've also tried to equilibrate a box of 216 water molecules with DFT and while the temperature and pressure seem to equilibrate quite easily, the constant of motion shows large jumps despite the potential and kinetic energy also being quite constant... If I'm not mistaken, the constant of motion in the NPT ensemble is the energy of the system + the energy of the thermostat + energy of barostat, as shown in this lecture from 2016, so this would suggest that the fluctuations are due to the thermostat or barostat, right? Like the TIP3P NPT run, I used a thermostat timecon of 100 fs and a barostat timecon of 1000 fs and after 21 125 fs (blue vertical line in figures), I restarted the calculation and increased the thermostat timecon to 500 fs. I accidentaly set RESTART_BAROSTAT and RESTART_BAROSTAT_THERMOSTAT to FALSE upon restarting, but I'm unsure if this might have lead to the large jumps in the constant of motion later at t > 28 000 fs... 

Same goes for the NPT_F calculation of my solvated TiO2 slab. I did a 10 ps simulation with thermostat timecon 50 and barostat timecon 100 and everything seems to equilibrate except for the constant of motion which seems to steadily rise. After 10 ps of simulation I restarted the calculation and again accidentaly messed up my input and set RESTART_BAROSTAT and RESTART_BAROSTAT_THERMOSTAT to FALSE but this time I didn't even change the timecons of the thermostat and barostat, but this time I immediately see a significant jump in the constant of motion after which it continues to increase... Any idea how I can fix this so that the constant of motion stays...well...constant?

 

Kind regards,

Léon

 

On Tuesday 23 January 2024 at 14:41:46 UTC+1 Léon Luntadila Lufungula wrote:

Dear Matthias,

 

My AIMD calculations are on a slab wherein I keep my bottom layer fixed, but apparently REGION MASSIVE is not compatible with constraints, so I'll have to use REGION GLOBAL. Either way, thanks a lot for the detailed answer! For a novice MD practicioner as myself this is already enough information to get started! I'll try the values you suggested during my production run while keeping an eye on the temperature and pressure during the simulation; I guess when I've chosen my TIMECON too large, this will result in the system drifting away from the equilibrium values.

 

As always, thanks a lot for all the help!

 

Kind regards,

Léon

On Tuesday 23 January 2024 at 14:21:05 UTC+1 Krack Matthias wrote:

Dear Léon

 

Using a MASSIVE CSVR thermostat with a small TIMECON of about 50 is appropriate and efficient for equilibration. Thereafter, the TIMECON for the thermostat can be increased to 500-1000 for the sampling (production) runs and also a GLOBAL instead of a MASSIVE thermostat can be used, if possible, to reduce the impact of the thermostat on the dynamics. You will find different opinions about the absolute and relative values for the TIMECONs of thermostat and barostat in the web or literature. I suggest to use a two to five times larger TIMECON for the barostat compared to the thermostat. You have, however, to keep in mind that the appropriate (best) parameters for baro- and thermostat depend on the type of system and the selected conditions eventually.

 

HTH

 

Matthias  

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Tuesday, 23 January 2024 at 11:47
To: cp2k <cp...@googlegroups.com>
Subject: Re: [CP2K:19801] Cons. Qty for different ensembles + Pressure fluctuation

Hi Matthias,

 

So I started an NpT directly from the inputs you gave me with a box containing 216 water molecules. After a 300ps equilibration (plots thus starrt at 0.3 ns), the pressure fluctuations are again in the 10kbar range which is okay as you pointed out, while the average pressure is around -14 bar, so tha'ts also within the 20 bar range you mentioned. The energy drift per atom is also below 1K so everything seems to work reasonably. 

 

I'm starting my AIMD simulations as we speak, so hopefully all goes well there too. I read that it is best to do a short equilibration with a strong coupling to the thermostat (e.g. CSVR, TIMECON 20 and REGION MASSIVE), before going to a weaker coupling (e.g. CSVR, TIMECON 200 and REGION GLOBAL), but is there also a rule of thumb w.r.t. the barostat TIMECON? I regularly see TIMECON 1000, but is it useful to set it to TIMECON 100 in this early equilibration?

 

Thanks for all the help already, it really helps me a lot given the fact that I don't have a senior researcher with a lot of hands-on CP2K experience around!

 

Kind regards,

Léon

 

--

You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to

cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/c498aa90-e014-4484-869e-a94877dd9813n%40googlegroups.com.

[Léon Luntadila Lufungula]

Dear Matthias,

1. The values plotted are indeed the conserved quantity that is printed in the 6th column of the .ener file.
2. Since there are still relatively large fluctuations in the cell volume, should I then continue the equilibration until the cell volume is constant? I can confirm that the computational cost of AIMD is quite large, but I would rather avoid cutting corners, because I don't want deterioration of the sampling quality as you point out.
3. If changing the thermostat/barostat parameters will change the value of the constant of motion, can you still use the constant of motion as a measure of system equilibration when equilibrating in two parts: strong early equilibration (small timecons: T50, B100) + weaker late equilibration (larger timecons: T100, B1000)?
4. Regarding EPS_DEFAULT, I did test the energy conservation of an NVE trajectory for several combinations of EPS_SCF/EPS_DEFAULT on a slab without solvent and 1.0E-7/1.0E-10 (first picture) gave similar results as  1.0E-7/1.0E-12 (second picture), but perhaps I should indeed reduce it further to ensure stable dynamics.
   
   
5. I agree that a noticeable energy drift will always occur if you just simulate long enough, but I wouldn't expect such a large drift (0.08 a.u. for the water box and 0.4 a.u. for the solvated slab if you disregard the jump when changing the thermostat/barostat) to already pop up when simulating for only 20 ps...
6. For the solvated slab, I was testing how thick I needed to make the water layer to correctly model the structure of liquid water so I simulated an NPT_F trajectory for a water layer of 10, 20, 30 and 40 Å. The 10 Å water layer is the one I showed in the previous message with the large energy drift, but the 20 Å water layer does seem to show an equilibrating constant of motion (see below). The input parameters are all the same and the structures were all solvated in a similar way by adding a vacuum layer to the slab with ASE and solvating the vacuum layer with Gromacs (gmx solvate). Any idea why the 20 Å trajectory does seem to equilibrate very easily, while the 10, 30 and 40 Å trajectories all show these weird jumps in the constant of motion?
   
   

As always, thanks a lot for all the help and my apologies for the barrage of questions, these calculations are really hard to setup correctly...

Kind regards,

Léon 

[Léon Luntadila Lufungula]

Dear Matthias,

I tried continuing the calculation with the 20 Å water layer, because this was the only one that looked like it was going towards a point where the constant of motion would also be equilibrated. However, upon restarting the calculation directly from the .restart file and only adjusting the maximum number of steps, so that everything is restarted (including the barostat), I immediately see larger fluctuations in the pressure and a resulting drift of the constant of motion... (see below) The only things I can think of that are inherently different between the two calculations is the number/type of nodes I'm running on and the SCF_GUESS. The 1st calculation was ran on 4 normal nodes and the wavefunction for the first step was restarted from a .wfn file of a geometry optimization which was performed before the MD run, while the restarted calculation was on 2 high-memory nodes and the wavefunction was started from a superposition of atomic wavefunctions as I read in the exercise on AIMD that storing the wavefunction can slow down your MD run, so this was not stored during the AIMD run. I'm currently running a calculation with the same node setup as the first calculation, but I was wondering if it could be possible that these changes are due to the fact that I am not starting from the .wfn file of the 1st calculation? 

I can't do these calculations in one go as they take longer than the maximum wall time on my HPC cluster, so if I can't properly restart them I can't equilibrate my system, let alone do a production run...

Kind regards,

Léon

[Krack Matthias]

Hi Léon

 

Restarting a run on a different number of nodes should not change the results markedly for a well-behaved system. I suggest to write the .wfn restart file after each MD step to ensure a smooth continuation of the run after a restart. As I wrote already, 10 ps are a by far too short time period for pressure equilibration. You might need at least 5-10 times longer time periods.

 

Best

 

Matthias

 

From: cp...@googlegroups.com <cp...@googlegroups.com> on behalf of Léon Luntadila Lufungula <Leon.luntad...@uantwerpen.be>
Date: Tuesday, 30 January 2024 at 14:39
To: cp2k <cp...@googlegroups.com>
Subject: Re: [CP2K:19853] Cons. Qty for different ensembles + Pressure fluctuation

Dear Matthias,

 

I tried continuing the calculation with the 20 Å water layer, because this was the only one that looked like it was going towards a point where the constant of motion would also be equilibrated. However, upon restarting the calculation directly from the .restart file and only adjusting the maximum number of steps, so that everything is restarted (including the barostat), I immediately see larger fluctuations in the pressure and a resulting drift of the constant of motion... (see below) The only things I can think of that are inherently different between the two calculations is the number/type of nodes I'm running on and the SCF_GUESS. The 1st calculation was ran on 4 normal nodes and the wavefunction for the first step was restarted from a .wfn file of a geometry optimization which was performed before the MD run, while the restarted calculation was on 2 high-memory nodes and the wavefunction was started from a superposition of atomic wavefunctions as I read in the exercise on AIMD that storing the wavefunction can slow down your MD run, so this was not stored during the AIMD run. I'm currently running a calculation with the same node setup as the first calculation, but I was wondering if it could be possible that these changes are due to the fact that I am not starting from the .wfn file of the 1st calculation? 

I can't do these calculations in one go as they take longer than the maximum wall time on my HPC cluster, so if I can't properly restart them I can't equilibrate my system, let alone do a production run...

 

Kind regards,

Léon

On Monday 29 January 2024 at 12:24:28 UTC+1 Léon Luntadila Lufungula wrote:

Dear Matthias,

1.   The values plotted are indeed the conserved quantity that is printed in the 6th column of the .ener file.

2.   Since there are still relatively large fluctuations in the cell volume, should I then continue the equilibration until the cell volume is constant? I can confirm that the computational cost of AIMD is quite large, but I would rather avoid cutting corners, because I don't want deterioration of the sampling quality as you point out.

3.   If changing the thermostat/barostat parameters will change the value of the constant of motion, can you still use the constant of motion as a measure of system equilibration when equilibrating in two parts: strong early equilibration (small timecons: T50, B100) + weaker late equilibration (larger timecons: T100, B1000)?

4.   Regarding EPS_DEFAULT, I did test the energy conservation of an NVE trajectory for several combinations of EPS_SCF/EPS_DEFAULT on a slab without solvent and 1.0E-7/1.0E-10 (first picture) gave similar results as  1.0E-7/1.0E-12 (second picture), but perhaps I should indeed reduce it further to ensure stable dynamics.

5.   I agree that a noticeable energy drift will always occur if you just simulate long enough, but I wouldn't expect such a large drift (0.08 a.u. for the water box and 0.4 a.u. for the solvated slab if you disregard the jump when changing the thermostat/barostat) to already pop up when simulating for only 20 ps...

6.   For the solvated slab, I was testing how thick I needed to make the water layer to correctly model the structure of liquid water so I simulated an NPT_F trajectory for a water layer of 10, 20, 30 and 40 Å. The 10 Å water layer is the one I showed in the previous message with the large energy drift, but the 20 Å water layer does seem to show an equilibrating constant of motion (see below). The input parameters are all the same and the structures were all solvated in a similar way by adding a vacuum layer to the slab with ASE and solvating the vacuum layer with Gromacs (gmx solvate). Any idea why the 20 Å trajectory does seem to equilibrate very easily, while the 10, 30 and 40 Å trajectories all show these weird jumps in the constant of motion?

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/28bfa788-4f09-498f-b760-615702976dfdn%40googlegroups.com.

[Léon Luntadila Lufungula]

Hi Matthias,

Thanks for the help yet again! 

1) Regarding restarting the 20 Å run: 

- I managed to get it to restart without the pressure fluctuations increasing markedly (see below). There seemed to be something wrong with the original restart file, but I can't look into it because I accidentaly overwrote it... Either way, this problem is solved as I have written a script that restarts my calculations without causing a sudden increase in the pressure fluctuations.

- For future cases, I'll be sure to write the .wfn file to ensure a smooth continuation of the run after a restart.

2) Regarding the drifting constant of motion for the 10 Å structure:

- Increasing EPS_DEFAULT to 10E-12 did not seem to solve the problem...

However, by reducing the time constants to 10(T) and 20(B) fs the large jumps in the constant of motion do seem to be gone (for the moment).

3) Pressure and temperature equilibration:

- If you say that 10 ps is too short for pressure equilibration, I'll be sure to continue equilibrating for a longer time to ensure proper equilibration. However, looking at the pressure plot, I would say that the pressure is more or less equilibrated, but I'm guessing you mean that the volume fluctuations should also reduce to mark full equilibration, right?

- I'm a bit puzzled about the ratio between the thermostat and barostat time constant you mentioned in a previous reply. You stated that the barostat TIMECON is always some multiple LARGER than the thermostat TIMECON, while at the same time you say that the pressure equilibrates less quickly than the temperature does. Wouldn't it then be better to set the barostat TIMECON to a lower value during early equilibration or is there a reason why the barostat TIMECON can't be lower than the thermostat TIMECON? Stating it in a more practical way: can I increase the thermostat TIMECON once the temperature has equilibrated while keeping a low barostat TIMECON to save some time on my equilibration?

I'm guessing these will be some of the last questions regarding this topic as it is now just a matter of simulating long enough for everything to equilibrate and then looking into how I can save some time on this equilibration for future runs (which I will figure out on my own and won't bother you with). Thanks for all of the help these last few weeks, it really helped me a lot in getting these MD simulations up and running and for this you have my heartfelt gratitude!

Kind regards,

Léon

[Léon Luntadila Lufungula]

Dear Matthias,

Thanks for all the help over the past few weeks, I really appreciate it!

For anyone reading this thread in the future, 

I figured out the reason for the jumps in the constant of motion! It wasn't due to EPS_DEFAULT, TIMESTEP or any of the other settings mentioned in this thread, but due to the fact I had not defined a reference cell in a calculation where the cell dimensions change. This resulted in the FFT grid not being constant throughout the calculation which seemed to cause these weird jumps in the constant of motion. Once a reference cell of 10-20% larrger than the calculation cell was defined the jumps disappeared for all of my troublesome calculations, as shown in the picture below for the 10 Å structure mentioned in my previous message. I have attached my input file for anyone requiring it as an example.

Kind regards,

Léon