Monopole Dipole and Quadrupole Corrections in MAPbI3

[reza_namakian1]

Dear VASP Team,

I would greatly appreciate your opinion regarding the "Monopole, Dipole, and Quadrupole Corrections" in MAPbI3, particularly considering that the MA cations in this perovskite may carry dipole moments.

I am not sure whether I have used the correct tags in the INCAR file to model this material properly. Specifically, I am concerned about the following correction settings, and whether they are necessary:

IVDW = 13
LDIPOL = T
IDIPOL = 4
DIPOL = 0.5 0.5 0.5

Additionally, I am unsure if the cell size is adequate for this purpose.

I have attached the INCAR, POTCAR, and POSCAR files.

Thank you in advance for your time and advice!

Reza.

[ferenc_karsai]

If you are unsure, please always do test calculations. For example look at the desired observable (energy, forces etc.) and see how it changes by enabling/disabling a given parameter or by varying the parameter. Ideally it should converge to a given value (meaning by increasing the parameter beyond a given value the desired observable is not changing anymore). Always go from low computational cost/accuracy parameters to higher ones, so that you minimize the computational effort.

We didn't need to use monopole, dipole etc. corrections in our machine learning calculations for energies, forces and stresses, but I don't know if they can be important for other quantities.

Your cell size seems to be ok (we used similar sizes for the training structures to train our machine learning force fields), but I can't tell for sure in your case because it depends on the purpose. If you are unsure, please make convergence tests where you start from small POSCARs and increase them.

[reza_namakian1]

Thanks a lot for the suggestions, Ferenc!

As you mentioned, I relaxed the ion/volume with and without corrections. I systematically toggled each keyword, IDIPOL and LDIPOL, on and off. Regardless of the settings, the energy values became extremely high and failed to converge. This confirms that these corrections are not necessary for this model. However, I’m now planning to develop a machine learning force field (MLFF) for a slab configuration of the 2D version of this perovskite (please see the attached PDF), and I have some concerns listed below:

1. VASP documentation advises to “Use any of these tags only after pre-converging the orbitals without the LDIPOL tag,” and mentions, “The center of charge should be set in the INCAR file (DIPOL = center of mass).” I’m unsure how to implement this during an MD-NPT run. Do you have any suggestions?

(Note: I turned off the DIPOL tag for the slab, which caused the convergence of electronic minimization to slow down significantly.)

2. In the attached PDF, I examined the convergence of the vacuum thickness using LDIPOL = T, IDIPOL = 3, and DIPOL = 0.5 0.5 0.5 via single point calculation. VASP raised warnings for thicknesses between 1 to 5 Å:

“The minimum charge density times the volume of the cell along the axis of the dipole correction is larger than 1E-1, which could indicate that your work function is not accurate as there is no field-free region in your cell. Please consider either increasing the size of your cell along the dipole correction (vacuum dimension) or perhaps increasing the precision of your calculation.”

But what worries me the most is that the periodic images appear to screen or interact with each other quite strongly up to 10 Å, and the trend only begins to saturate around 50 Å. Do I really need to have that much vacuum? Also, do you think a cutoff radius of 5 or 6 Å would be sufficient based on this trend? I’m concerned that going beyond 7 Å might be challenging for the MLFF.

[ferenc_karsai]

I've talked to a colleague (Jonathan Lahnsteiner) who worked a lot with perovskites. He said he tried the dipol corrections for MAPbI3 and it never converged. Also the dipol corrections are not learned during on-the-fly MLFF. So don't use them.

We think the vacuum thickness plot you showed is probably ok for 10-20 Angstrom and would look even better if you would change the total energy to energy per atom.

For the parameters, there is a publication by Menno Bokdam who looked at similar things as you:
https://pubs.acs.org/doi/full/10.1021/acs.jpcc.3c01634
The calculation was done with the VASP MLFF so you could use the publication as a guideline for your parameters.

He used 8 and 5 Angstrom for the radial and angular cut-off, which are the default. So no need to change anything here.
What I saw is that he used an angular cutoff of 4. However, we changed the default to 3 since we realized that almost always an angular cutoff of 3 is sufficient. So please test this for your case.

[reza_namakian1]

Ferenc,

I wanted to provide a quick update on the dipole correction. I turned off the dipole tags and repeated the vacuum thickness convergence analysis. I noticed that not using the dipole correction yields even slightly lower energies for the entire range of thicknesses I tested before. This makes me more certain that the dipole correction might not improve the current model I'm working on.

Additionally, I’ve started experimenting with the ML tags on a simple system involving an MA cation as a toy model (please see the attached zipped file for more details and files). I plan to extend this MLFF to include PbI octahedra and, eventually, the larger cations like Butylammonium (BA).

A few observations and questions:

1) I attempted to rename the POSCAR file header to “MA Cation” based on the "Mind" box provided in the documentation for "ML_IWEIGHT = 3". However, I noticed the following in the log file:
*Use also structure names for distinction into groups for weighting: F ML_LUSE_NAMES*

2) Despite explicitly setting ML_LSPARSDES to `.TRUE.` in the INCAR, the log file indicates:
*Enable sparsification of angular descriptors: F (!) ML_LSPARSDES*

3) The log file also notes:
*Specifies whether self-interactions are subtracted or not: T (?) ML_LSIC*

4) Lastly, following the documentation, which suggests creating distinct groups for different species (e.g., O1 and O2 for oxygen in bulk vs. surface molecules), I’m considering distinguishing MA and BA species in the POSCAR files. For MA, I intend to label them as N1, C1, and H1, and for BA, as N2, C2, and H2. Am I interpreting this correctly, or have I misunderstood this approach?

Thank you in advance for your advice!

[ferenc_karsai]

So to the four points you have observed:

1) The code by default divides the substructures into groups containing the exact same types and atoms per type. If you want to further subdivide the structures by giving names you have to set ML_LUSE_NAMES=.TRUE.. The description of ML_IWEIGHT was slightly wrong and I changed the description and added a wiki page for ML_LUSE_NAMES. ML_LUSE_NAMES is only important if you have structures with significantly different energies, e.g. structures at vastly different pressures but still the same atom types and number of atoms per type. I don't think you have this in your example, so you don't have to bother with this tag.

2) ML_LSPARSDES is only implemented for the ML_MODE=refit and refitbayesian. You tried to use it during ML_MODE=train. In VASP.6.4.3 the code will automatically overwrite this tag. I added a description on the wikipage (https://www.vasp.at/wiki/index.php/ML_LSPARSDES) that this mode is only available for ML_MODE=refit and refitbayesian and in the next release version the code will give a graceful exit when trying to use it with ML_MODE=train.
I think you should not use ML_LSPARSDES but rather use ML_DESC_TYPE=1. This will make your force field 10-20% less accurate but around 3.5 times faster, which in my opinion is in most cases a good tradeoff. But beware the same is true for ML_DESC_TYPE, that you have to use it only after your ML_MODE=train calculations are finished and you retrain with ML_MODE=refit.

3) This tag is only there for debugging purposes for developers and the user should not bother with this tag. It is also on purpose not described.

4) Yes you are interpreting correctly. However, I asked again my colleague who has huge experience with MLFFs for perovskites and he said they did similar calculations but didn't have to distinguish between different groups of Oxygens and Nitrogens. This is important because the machine learning force field in the production runs scales quadratically with the number of species (linearly for some parts ith ML_DESC_TYPE=1). So if you plan to do massive production calculations, it is maybe better to train without distinction of different Oxygens, etc. and test if this force field is good enough. If it is not good enough then modify your ML_AB file to distinct types (this will require some scripting), select the local reference configurations with the new species (ML_MODE=select) and then refit (ML_MODE=refit).

[reza_namakian1]

Many thanks, Frence! I learned a lot from this post.

I have just one last question. I noticed the very useful feature of "Merging different ML_AB files." This is quite helpful for cases where I want for example to experiment with replacing cations or other molecules and merge them with separately already trained structures. I was wondering if you guys have any script (or perhaps a module in py4vasp) to automate this merging process, as explained in the documentation, to reduce human interference and minimize errors during this process.

[ferenc_karsai]

Sorry, but we have no scripts that are shareable (but we are planning modules that do that). But if possible try to avoid merging separate calculations, since the ML_MODE=select can become quite time consuming and if you have to merge separate files you have to be very careful to not make a mistake.
I would advise to first train on bulk, then continue training the ML_AB file from that calculation on the surface than use the ML_AB from there for the next training and so on.

If you still insist on merging files be very careful to set entries like "The maximum number of atoms per system", "The maximum number of atoms per atom type", etc. correctly among all systems. This is something where users sometimes make mistakes.
If you reselect with ML_MODE=select you don't have to set the number of local reference configurations correctly. Just set 1 for each and "1 1" for "Basis set for X".