MD calculation of NPT_F ensemble including baro- and thermostats at various T

[Lukas C]

Hi everyone!

I am trying to see signs of melting for a given SiO2 bulk. In order to achieve this, I ran MD calculations for a NPT_F (constant pressure and temperature, flexible cell) ensemble which is connected to a thermal bath (adaptive langevin thermostat) and to a barostat which are utilized to keep T and p at constant values for each simulation respectively. The temperature was set to T=600K,900K,1200K,1500K,1800K,2100K and 3000K.

However, i am experiencing multiple problems.
1. The calculations abort after 20 steps (every time) with the error:

 *******************************************************************************
 *   ___                                                                       *
 *  /   \                                                                      *
 * [ABORT]                                                                     *
 *  \___/                  invalid value for enumeration:   104                *
 *    |                                                                        *
 *  O/|                                                                        *
 * /| |                                                                        *
 * / \                                     input/input_enumeration_types.F:188 *

 *******************************************************************************

2. by comparing the output files, i found out that the simulation seems to trigger the same reactions (altough the temperature is very different). cp2k did recognice the temperature (thus, the kinetic energy is changed) but the curves that represent the change of kinetic energy resemble each other perfectly. to my understanding, this means that exactly the same reactions were triggered, which should not be the case for different temperatures. Each simulation was started with a .wfn file as an initial guess, which was taken from an single-point ENERGY calculation that was conducted beforehand. Thus the potential energy is the same for each simulation at t=0. The plot is attached together with the input file. I already tried to use various optimizers, preconditioners and so on. am i doing something seriously wrong? Thank you, your help is appreciated!

Regards,
Lukas

[Marcella Iannuzzi]

Dear Lukas, 

Did you  first  equilibrate the system at room temperature?  

Why did you choose the Langevin thermostat? How does the cell fluctuate? Which is the energy of the barostat? The stress tensor?

Is the total energy conserved? 

Without more information it is hard to guess what the problem might be. 

If you have an initially equilibrated sample (stress tensor around zero, conserved total energy, etc. ), try annealing it progressively to higher temperatures.

The model looks quite small and it will be anyway difficult to get a disordered structure due to the strong spatial correlations. 

Kind regards

Marcella 

[Lukas C]

Hi Marcella! Hi everyone!

EDIT: i cannot post with files attached. i will try to do it seperately.

Many thanks for your quick reply! Let me try to answer your questions.

1. The structure that i am using as an input was made by means of the melt and quench technique with ReaxFF force field MD in LAMMPS. After that it was cell-optimized in CP2K at 0K. After relaxation, the entries S_ij of the stress tensor satisfied S_ij<10MPa.
2. The langevin thermostat was only used because it did not seem wrong to use it. I see that the barostat has also a section called thermostat. should i just use this one?

3. The total energy is not concerved since I am using an NPT (Gibbs) ensemble.
4 .The cell fluctuations are given in the attached file (exemplary for T=3000K).
5 .I did not print the barostat-energy until now but I will send this information as soon as possible, most probably by tomorrow.

Concerning your advice to try annealing:
I did already try to do that last week. I simply used this as an &MD input:
 &MD
    ENSEMBLE NVE        
    INITIALIZATION_METHOD DEFAULT
    MAX_STEPS  1000
    STEPS      1000
    TIMESTEP 0.5        
    TEMPERATURE 300
    ANNEALING 1.001
 &END MD

and this seemed to work in principle but it ran into convergence problems... please see the attached output file (I removed stuff from the middle of the file to shorten it a bit).

Thank you! Your help is greatly appreciated.

Regards,
Lukas

[Lukas C]

somehow sending failes if something is attached. so here is at least the cell file directly:

#   Step   Time [fs]       Ax [Angstrom]       Ay [Angstrom]       Az [Angstrom]       Bx [Angstrom]       By [Angstrom]       Bz [Angstrom]       Cx [Angstrom]       Cy [Angstrom]       Cz [Angstrom]      Volume [Angstrom^3]
       0       0.000       14.6984880000        0.0000000000        0.0000000000        0.3089790000       14.7861450000        0.0000000000        0.1985100000        0.5067420000       15.2248250000          3308.8717336268
       1       0.500       14.7497158555        0.0537881155       -0.0746164496        0.3641647584       14.7132968826        0.0232584917        0.1237678600        0.5304966664       15.3174452448          3323.7848176730
       2       1.000       14.8203073638        0.0909489367       -0.1278713294        0.4030562136       14.6928503265        0.0408367314        0.0708259428        0.5495397021       15.4169421853          3356.2858676037
       3       1.500       14.8815117820        0.0998885160       -0.1403322084        0.4134637396       14.7157181601        0.0472078290        0.0589874203        0.5572786312       15.4804000421          3389.1447391791
       4       2.000       14.9060983270        0.0799887119       -0.1042040876        0.3942208231       14.7646399509        0.0435861959        0.0958924887        0.5542243288       15.4751254001          3405.0917706236
       5       2.500       14.8812900116        0.0403580450       -0.0270083285        0.3541585911       14.8191072238        0.0314274788        0.1739252767        0.5414117396       15.3970929978          3395.0727563140
       6       3.000       14.8138745269       -0.0061503345        0.0670720618        0.3062662026       14.8699609176        0.0123059514        0.2686233223        0.5208331067       15.2763705366          3364.7820027850
       7       3.500       14.7337595627       -0.0450625286        0.1468154298        0.2657252633       14.9200057624       -0.0071032058        0.3485556378        0.5003642963       15.1596087586          3331.9931193541
       8       4.000       14.6773702385       -0.0638454147        0.1887667953        0.2459403984       14.9709452176       -0.0188316933        0.3903565790        0.4889183411       15.0949860284          3316.1754306774
       9       4.500       14.6682754715       -0.0582819761        0.1863771395        0.2516633294       15.0203060127       -0.0174964551        0.3876777063        0.4920873408       15.1088612241          3328.1004091840
      10       5.000       14.7075619173       -0.0305655318        0.1424180549        0.2807217635       15.0550067870       -0.0027367320        0.3433762601        0.5096603743       15.1922818057          3363.3472150915
      11       5.500       14.7735852436        0.0133400318        0.0704813688        0.3267183782       15.0495825368        0.0187897572        0.2710381271        0.5336189236       15.3103513690          3403.5567335316
      12       6.000       14.8367420057        0.0613415460       -0.0063946166        0.3768712692       14.9871874829        0.0398887289        0.1937403565        0.5553211906       15.4154820085          3427.1352009275
      13       6.500       14.8755539262        0.0996359595       -0.0661684631        0.4168405581       14.8769757073        0.0531995258        0.1335493133        0.5668240425       15.4661719421          3421.7398071290
      14       7.000       14.8884917853        0.1176999613       -0.0891680748        0.4360179505       14.7536209773        0.0543571012        0.1103968944        0.5641980157       15.4523027199          3393.1143082129
      15       7.500       14.8866971314        0.1115060801       -0.0659352382        0.4305733264       14.6614395092        0.0430280764        0.1340423708        0.5482800246       15.3947014767          3359.0781832207
      16       8.000       14.8826831602        0.0862109392       -0.0049867059        0.4057660296       14.6424113987        0.0210144166        0.1959968464        0.5226053131       15.3273176472          3339.4180626402
      17       8.500       14.8847638718        0.0521850999        0.0736710963        0.3720557197       14.7101532776       -0.0063911550        0.2759903128        0.4937283449       15.2850165172          3346.2283225409
      18       9.000       14.8856128696        0.0215935376        0.1442514111        0.3414763728       14.8415897429       -0.0312917259        0.3475748927        0.4698799980       15.2806271438          3375.2752408601
      19       9.500       14.8672397840        0.0000282921        0.1845898569        0.3193906853       14.9884536842       -0.0452332231        0.3876174048        0.4587965656       15.3087700030          3410.6224049567
      20      10.000       14.8187679059       -0.0088967582        0.1846394360        0.3093192078       15.0934343643       -0.0426562123        0.3856188882        0.4645482911       15.3498769842          3432.5350392072

[Marcella Iannuzzi]

Dear Lukas, 

It is recommendable to monitor the behaviour of the system at conditions closer to equilibrium. 

Jumping from 0K to 3000K seems quite ambitious. 

The conserved quantity has to be conserved also with the NPT ensemble, where it is not simply the sum of pot. and kin. energy, but it includes all the due terms of the extended ensemble. If it is not conserved there is a problem.  

I would be conservative with the choice of the thermostat, Nose-Hoover and CSVR are for sure reliable. 

Adding a thermostat to the barostat helps in controlling the cell fluctuations, that indeed are quite wild according to the posted example. 

Anyway, monitoring the stress is also useful to check on the sanity of the system. 

For the annealing, I would rather suggest a step-wise procedure, i.e. running a series of run at increasing temperature  (\Delta T ~ 200 K ), starting always from the previous one. 

Kind regards

Marcella

[Lukas C]

Dear Marcella!

Thank you so much for your help! I understand a lot more now. I will try to incorporate your advices and keep you updated about any progress.

Best regards,
Lukas

[Lukas C]

Hi Marcella! Hi everyone!

I tried to run the annealing as suggested. So I turned off the thermostats and choose an NPE_F ensemble. The temperature is changed very gently: starting from 2K, with a rescaling factor for the velocities of 1.01, it should reach 3.127K after 50 steps. Unfortunately, I get unexplained MPI errors at the end of the first SCF step which apart from that converges nicely and does not show any signs of malfunction. I do get a WARNING before the SCF cycle starts, namely when calculating the PDOS at iteration step 0. All the input and output files (except for the .wfn file which is too large) are attached. I tried running it with different preconditioners but this did not change anything. I even ran the simulation on another cluster where it got stuck at the same point.
Any thoughts on this? Thank you!

Kind regards,
Lukas

[Lukas C]

Hi!

I just found out that it works with an NVE ensemble. but this seems to be a bad choice in my understanding.

Best regards,
Lukas