Liquid Si - MLFF

From VASP Wiki

Task

Generating a machine learning force field for liquid Si. For this tutorial, we expect that the user is already familiar with running conventional ab initio molecular dynamic calculations.

Input

POSCAR

In this example we start from a 64 atom super cell of diamond-fcc Si (the same as in Liquid Si - Standard MD):

Si cubic diamond 2x2x2 super cell of conventional cell
     5.43090000000000
    2.00000000   0.00000000   0.00000000
    0.00000000   2.00000000   0.00000000
    0.00000000   0.00000000   2.00000000
   Si
   64
Direct
   0.00000000   0.00000000   0.00000000
   0.50000000   0.00000000   0.00000000
   0.00000000   0.50000000   0.00000000
   0.50000000   0.50000000   0.00000000
   0.00000000   0.00000000   0.50000000
   0.50000000   0.00000000   0.50000000
   0.00000000   0.50000000   0.50000000
   0.50000000   0.50000000   0.50000000
   0.37500000   0.12500000   0.37500000
   0.87500000   0.12500000   0.37500000
   0.37500000   0.62500000   0.37500000
   0.87500000   0.62500000   0.37500000
   0.37500000   0.12500000   0.87500000
   0.87500000   0.12500000   0.87500000
   0.37500000   0.62500000   0.87500000
   0.87500000   0.62500000   0.87500000
   0.00000000   0.25000000   0.25000000
   0.50000000   0.25000000   0.25000000
   0.00000000   0.75000000   0.25000000
   0.50000000   0.75000000   0.25000000
   0.00000000   0.25000000   0.75000000
   0.50000000   0.25000000   0.75000000
   0.00000000   0.75000000   0.75000000
   0.50000000   0.75000000   0.75000000
   0.37500000   0.37500000   0.12500000
   0.87500000   0.37500000   0.12500000
   0.37500000   0.87500000   0.12500000
   0.87500000   0.87500000   0.12500000
   0.37500000   0.37500000   0.62500000
   0.87500000   0.37500000   0.62500000
   0.37500000   0.87500000   0.62500000
   0.87500000   0.87500000   0.62500000
   0.25000000   0.00000000   0.25000000
   0.75000000   0.00000000   0.25000000
   0.25000000   0.50000000   0.25000000
   0.75000000   0.50000000   0.25000000
   0.25000000   0.00000000   0.75000000
   0.75000000   0.00000000   0.75000000 
   0.25000000   0.50000000   0.75000000
   0.75000000   0.50000000   0.75000000
   0.12500000   0.12500000   0.12500000
   0.62500000   0.12500000   0.12500000
   0.12500000   0.62500000   0.12500000
   0.62500000   0.62500000   0.12500000
   0.12500000   0.12500000   0.62500000
   0.62500000   0.12500000   0.62500000
   0.12500000   0.62500000   0.62500000
   0.62500000   0.62500000   0.62500000
   0.25000000   0.25000000   0.00000000
   0.75000000   0.25000000   0.00000000
   0.25000000   0.75000000   0.00000000
   0.75000000   0.75000000   0.00000000
   0.25000000   0.25000000   0.50000000
   0.75000000   0.25000000   0.50000000
   0.25000000   0.75000000   0.50000000
   0.75000000   0.75000000   0.50000000
   0.12500000   0.37500000   0.37500000
   0.62500000   0.37500000   0.37500000
   0.12500000   0.87500000   0.37500000
   0.62500000   0.87500000   0.37500000
   0.12500000   0.37500000   0.87500000
   0.62500000   0.37500000   0.87500000
   0.12500000   0.87500000   0.87500000
   0.62500000   0.87500000   0.87500000

KPOINTS

We will start with a single k point in this example:

K-Points
 0
Gamma
 1  1  1
 0  0  0

INCAR

#Basic parameters
ISMEAR = 0
SIGMA = 0.1
LREAL = Auto
ISYM = -1
NELM = 100
EDIFF = 1E-4
LWAVE = .FALSE.
LCHARG = .FALSE.

#Parallelization of ab initio calculations
NCORE = 2

#MD
IBRION = 0
MDALGO = 2
ISIF = 2
SMASS = 1.0
TEBEG = 2000
NSW = 10000
POTIM = 3.0
RANDOM_SEED =          88951986                0                0

#Machine learning paramters
ML_FF_LMLFF = .TRUE.
ML_FF_ISTART = 0
ML_FF_LMLFF = .TRUE.
switches on the machine learning of the force field
ML_FF_ISTART = 0
corresponds to learning from scratch

Calculation

Creating the liquid structure

Because we don't have a structure of liquid silicon readily available, we first create that structure by starting from a super cell of crystalline silicon with 64 atoms. The temperature is set to 2000 K so that the crystal melts rapidly in the MD run. To improve the simulation speed drastically, we utilize the on-the-fly machine learning. Most of the ab initio steps will be replaced by very fast force-field ones. Within 10000 steps equivalent to 30 ps, we have obtained a good starting position for the subsequent simulations in the CONTCAR file. You can copy the input files or download them.

After running the calculation, we obtained a force field, but its initial trajectory might be tainted but the unreasonable starting position. Nevertheless, it is instructive to inspect the output to understand how to assess the accuracy of a force field, before refining it in subsequent calculations. The main output files for the machine learning are

ML_ABN
contains the ab initio structure datasets used for the learning. It will be needed for continuation runs as ML_AB.
ML_FFN
contains the regression results (weights, parameters, etc.). It will be needed for continuation runs as ML_FF.
ML_LOGFILE
logging the proceedings of the machine learning. This file consists of keywords that are nicely "grepable." The keywords are explained in the in the beginning of the file and upon "grepping". The status of each MD step is given by the keyword "STATUS". Please invoke the following command:
grep STATUS ML_LOGFILE

The output should look similar to the following:

# STATUS ###############################################################
# STATUS This line describes the overall status of each step.
# STATUS 
# STATUS nstep ..... MD time step or input structure counter
# STATUS state ..... One-word description of step action
# STATUS             - "accurate"  (1) : Errors are low, force field is used
# STATUS             - "threshold" (2) : Errors exceeded threshold, structure is sampled from ab initio
# STATUS             - "learning"  (3) : Stored configurations are used for training force field
# STATUS             - "critical"  (4) : Errors are high, ab initio sampling and learning is enforced
# STATUS             - "predict"   (5) : Force field is used in prediction mode only, no error checking
# STATUS is ........ Integer representation of above one-word description (integer in parenthesis)
# STATUS doabin .... Perform ab initio calculation (T/F)
# STATUS iff ....... Force field available (T/F, False after startup hints to possible convergence problems)
# STATUS nsample ... Number of steps since last reference structure collection (sample = T)
# STATUS ngenff .... Number of steps since last force field generation (genff = T)
# STATUS ###############################################################
# STATUS            nstep     state is doabin    iff   nsample    ngenff
# STATUS                2         3  4      5      6         7         8
# STATUS ###############################################################
STATUS                  0 threshold  2      T      F         0         0
STATUS                  1 critical   4      T      F         0         1
STATUS                  2 critical   4      T      F         0         2
STATUS                  3 critical   4      T      T         0         1
STATUS                  4 critical   4      T      T         0         1
STATUS                  5 critical   4      T      T         0         1
     .                  .        .   .      .      .         .         .
     .                  .        .   .      .      .         .         .
     .                  .        .   .      .      .         .         .
STATUS               9997 accurate   1      F      T       945       996
STATUS               9998 accurate   1      F      T       946       997
STATUS               9999 accurate   1      F      T       947       998
STATUS              10000 learning   3      T      T       948       999

Another important keyword is "ERR". For this instance we should type the following command:

grep ERR ML_LOGFILE

The output should look like the following:

# ERR ######################################################################
# ERR This line contains the RMSEs of the predictions with respect to ab initio results for the training data.
# ERR 
# ERR nstep ......... MD time step or input structure counter
# ERR rmse_energy ... RMSE of energies
# ERR rmse_force .... RMSE of forces
# ERR rmse_stress ... RMSE of stress
# ERR ######################################################################
# ERR               nstep      rmse_energy       rmse_force      rmse_stress
# ERR                   2                3                4                5
# ERR ######################################################################
ERR                     0   0.00000000E+00   0.00000000E+00   0.00000000E+00
ERR                     1   0.00000000E+00   0.00000000E+00   0.00000000E+00
ERR                     2   0.00000000E+00   0.00000000E+00   0.00000000E+00
ERR                     3   2.84605192E-05   9.82351889E-03   2.40003743E-02
ERR                     4   1.83193349E-05   1.06700600E-02   5.37606479E-02
ERR                     5   4.12132223E-05   1.34123085E-02   1.01588957E-01
ERR                     6   9.51627413E-05   1.90335214E-02   1.31959103E-01
.                       .                .                .                .
.                       .                .                .                .
.                       .                .                .                .
ERR                  9042   1.07159240E-02   2.41283323E-01   4.95695745E+00
ERR                  9052   1.07159240E-02   2.41283323E-01   4.95695745E+00
ERR                 10000   1.07159240E-02   2.41283323E-01   4.95695745E+00

This tag lists the errors on the energy, forces and stress of the force field compared to the ab initio results on the available training data. The second column shows the MD step. We see that the entry is not output at every MD step. The errors only change if the force field is updated, hence when an ab initio calculation is executed (it should correlate with the doabin column of the STATUS keyword). The other three columns show the errors on the energy (eV/atom), forces (ev/Angstrom) and stress (kB).

Structral properties of the force field

To examine the accuracy of structural properties, we compare the deviations between a 3 ps molecular dynamics run using the force field and a full ab initio calculation. For a meaningful comparison, it is best to start from the same initial structure. We will use the liquid structure, we obtained in the previous step and back it up

cp CONTCAR POSCAR.T2000_relaxed

Now, we proceed with the force field calculation and set up the required files

cp POSCAR.T2000_relaxed POSCAR
cp ML_FFN ML_FF

To run a shorter simulation using only the force field, we change the following INCAR tags to

ML_ISTART = 2
NSW = 1000

After the calculation finished, we backup the history of the atomic positions

cp XDATCAR XDATCAR.MLFF_3ps

To analyze the pair correlation function, we use the PERL script pair_correlation_function.pl

Click to show/hide pair_correlation_function.pl

and process the previously saved XDATCAR files

perl pair_correlation_function.pl XDATCAR.MLFF_3ps > pair_MLFF_3ps.dat

To save time the pair correlation function for 1000 ab initio MD steps is precalculated in the file pair_AI_3ps.dat.

The interested user can of course calculate the results of the ab initio MD by rerunning the above steps while switching off machine learning via

ML_LMLFF = .FALSE.

We can compare the pair correlation functions, e.g. with gnuplot using the following command

gnuplot -e "set terminal jpeg; set xlabel 'r(Ang)'; set ylabel 'PCF'; set style data lines; plot 'pair_MLFF_3ps.dat', 'pair_AI_3ps.dat' " > PC_MLFF_vs_AI_3ps.jpg

The pair correlation functions obtained that way should look similar to this figure

We see that pair correlation is quite well reproduced although the error in the force of 0.247664 eV/ shown above is a little bit too large. This error is maybe too large for accurate production calculations (usually an accuracy of approximately 0.1 eV/ is targeted), but since the pair correlation function is well reproduced it is perfectly fine to use this on-the-fly force field in the time-consuming melting of the crystal.

Obtaining a more accurate force field

Including the melting phase in the force field may impact the accuracy of the force field. To improve it is usually advisable to learn on the pure structures, which in our case this means to use the CONTCAR file obtained after the melting. Furthermore, the force field was learned using only a single k point so that we can improve the accuracy of the reference data by including more k points. In most calculations, it is important to conduct accurate ab initio calculations since bad reference data can limit the accuracy or even inhibit the learning of a force field.

We restart from the liquid structure obtained before

cp POSCAR.T2000_relaxed POSCAR

and change the following INCAR tags

ALGO = Normal
ML_LMLFF = .TRUE.
ML_ISTART = 0
NSW = 10000

If you run have resources to run in parallel, you can reduce the computation time further by adding k point parallelization with the KPAR tag. We use a denser k-point mesh in the KPOINTS file

2x2x2 Gamma-centered mesh
0 0 0
Gamma
 2 2 2
 0 0 0

We will learn a new force field with 1000 MD steps (each of 3 fs). Keep in mind to run the calculation using the standard version of VASP (usually vasp_std). After running the calculation, we examine the error "ERR" in the ML_LOGFILE by typing:

grep "ERR" ML_LOGFILE

where the last entries should be close to

ERR                   886   5.98467749E-03   1.48190308E-01   2.38264786E+00
ERR                   908   5.98467749E-03   1.48190308E-01   2.38264786E+00
ERR                   925   5.98467749E-03   1.48190308E-01   2.38264786E+00
ERR                   959   5.98467749E-03   1.48190308E-01   2.38264786E+00
ERR                   980   5.98467749E-03   1.48190308E-01   2.38264786E+00
ERR                   990   5.99559653E-03   1.50261779E-01   2.40349561E+00
ERR                  1000   5.99559653E-03   1.50261779E-01   2.40349561E+00

We immediately see that the errors for the forces are significantly lower than in the previous calculation with only one k point. This is due to the less noisy ab initio data which is easier to learn.

To understand how the force field is learned, we inspect the ML_ABN file containing the training data. In the beginning of this file, you will find information about the number of reference structures

 1.0 Version
**************************************************
     The number of configurations
--------------------------------------------------
         48

and the number of local reference configurations (size of the basis set)

**************************************************
     The numbers of basis sets per atom type
--------------------------------------------------
       382

We will monitor how much training data is added to a structure after each learning cycle and what impact this has on the accuracy of the force field. First we will perform a continuation run keeping the Bayesian threshold fixed

cp ML_ABN ML_AB
cp CONTCAR POSCAR

and set the following INCAR tags

ML_ISTART = 1
ML_CTIFOR = x

where x is the threshold of the Bayesian error obtained from the ML_LOGFILE by typing

grep BEEF ML_LOGFILE

and using the value for the threshold (6th column) in the last entry. By setting a good estimate for the threshold several ab initio steps are skipped in the first few steps, which would be required to determine the threshold automatically. Mind: When continuation runs are performed for different crystal structures, the last previous threshold for the Bayesian error might be too large leading to premature skipping of ab initio steps. This is particular relevant when studying the liquid phase, first, and applying its threshold to the solid phase. In that case it is safer to use the default value for ML_CTIFOR = .

After running the calculation, we inspect the last instance of the errors in the ML_LOGFILE, which should look similar to

====================================================================================================
     Information on error estimations
----------------------------------------------------------------------------------------------------
              Error in energy (eV atom^-1):     0.006252
             Error in force  (eV Angst^-1):     0.168549
                      Error in stress (kB):     2.606332
               Bayesian error (eV Angst-1):     0.037569    0.101357
                       Spilling factor (-):     0.001114    0.020000
====================================================================================================

We see that the accuracy has not changed significantly and the threshold for the Bayesian error is set to the one specified in the INCAR file. Looking at the actual errors we observe that the threshold was exceeded only a few times. You can find all these instances with the command

grep "Bayesian error (eV Angst-1):" ML_LOGFILE | awk '{if ($5 > $6) { print $5 }}'

Consequently only very few structures were added to the ab initio data in the ML_ABNCAR file. The entry for the reference structures and basis sets should look similar to

**************************************************
     The number of configurations
--------------------------------------------------
         65
...
**************************************************
     The numbers of basis sets per atom type
--------------------------------------------------
       407

From this observation, we conclude either that the force field is already accurate enough or that the threshold is set too high. In the latter case, the machine would not conduct ab initio calculations because the force field is judged accurate enough according to the threshold.

We compare this behavior to the case of using the default value for ML_FF_CTIFOR. We continue from the existing structure and force field

cp CONTCAR POSCAR
cp ML_ABNCAR ML_ABCAR

Remove the ML_FF_CTIFOR tag from the INCAR file.

The default value for ML_FF_CTIFOR is practically zero, so we expect more structures to be added. We inspect the ML_ABNCAR to investigate how the number of structures and the basis set size changed

 **************************************************
     The number of configurations
--------------------------------------------------
        160
...
**************************************************
     The numbers of basis sets per atom type
--------------------------------------------------
       726

Even though the number of structure datasets and the basis set size increased, the accuracy of the force field did not improve significantly

====================================================================================================
     Information on error estimations
----------------------------------------------------------------------------------------------------
              Error in energy (eV atom^-1):     0.005950
             Error in force  (eV Angst^-1):     0.162064
                      Error in stress (kB):     2.628567
               Bayesian error (eV Angst-1):     0.029723    0.061348
                       Spilling factor (-):     0.000406    0.020000
====================================================================================================

Ideally, one should continue learning until no structures need to be added to the training data and basis set anymore. In practice, the prediction of the Bayesian error exhibits numerical inaccuracies so that an ab initio calculation is conducted from time to time even if the force field is accurate enough. Hence, in this tutorial, we stop at this point because further addition of training data does not improve the accuracy on the training structures compared to the ab initio data.

Download

MLFF_Liquid_Si.tgz