My Community

Dear VASP Community,

I am currently working on a complex system involving a self-assembled monolayer of BPDCA (C14H10O4) adsorbed on a CoO(100) surface, which is modeled with 3 layers of CoO(100) (totaling 108 atoms).

After running the geometry optimization, the calculated adsorption energy for this system appears to be unreasonably high, specifically -25.770732 eV. To ensure that the issue isn't related to the individual components of the system, I revisited the parameters used in the substrate calculations. I conducted tests by adjusting the SIGMA parameter (for Gaussian smearing) from its original value of 0.1 eV to 0.01 eV. Notably, the ionic optimization of the substrate using SIGMA = 0.01 eV resulted in an unexpectedly low energy for the CoO (100) substrate, around -641.07897822 eV, which will significantly alters the adsorption energy by approximately 14.45 eV. However, to test the reliability of the SIGMA parameter and to make sure that I am looking at the true ground state energy for the substrate. I ran an additional single point calculations just for the substrate.

Here are the results from my scf for different SIGMA:

SIGMA = 0.01 eV: -626.13568056 eV
SIGMA = 0.02 eV: -627.97482516 eV
SIGMA = 0.03 eV: -626.15113909 eV

Additionally, when I altered the MAGMOM parameter from MAGMOM = 27*4.0 27*-4.0 54*0.0 to MAGMOM = 27*5.0 27*-5.0 54*0.0, the SCF total energy further decreased to -646.03208699 eV. These tests suggest that the energy -641.07897822 eV cannot be the true ground state of the CoO(100) substrate, indicating that the system's energy landscape is highly sensitive to even minor changes in parameters. Consequently, determining the correct electronic and magnetic convergence is crucial.

Moreover, in additional tests aimed at ensuring proper convergence, I reduced the SCF convergence criterion to 10^-8 eV and experimented with charge density and magnetic moment mixing parameters. The default settings for AMIX, AMIX_MAG, BMIX, and BMIX_MAG resulted in a substrate energy of -642.30933566 eV. However, following the recommended settings from the VASP Wiki (AMIX = 0.2, BMIX = 0.0001, AMIX_MAG = 0.8, BMIX_MAG = 0.0001), the SCF energy was -627.97483064 eV. Further reducing AMIX to 0.1 and AMIX_MAG to 0.4 gave me an electronic and magnetic ground state identical to the one found with the MAGMOM = 27*5.0 27-5.0 54*0.0 setting, with a TOTEN of -646.03206512 eV. Interestingly, altering the SIGMA value in these new tests did not significantly change the results. I re-optimized my CoO(100) with the AMIX and AMIX_MAG that yields -646.03206512 eV and the new TOTEN I get with these parameters is -646.56425475 eV and the scf energy for this job is -646.56425554 eV. Does it mean that AMIX=0.1 and AMIX_MAG=0.4 with nearly zero BMIX and BMIX_MAG gives me the lowest poosible energy state of CoO(100) substate ? And I could consider this as the correct magnetic ground state.

I have attached the zip of the original input files for both the BPDCA on CoO(100) system and the CoO(100) substrate optimization and the MAIN output files. Additionally, I have attached the OUTCAR for the scf job of AMIX=0.1 and AMIX_MAG=0.4 job with name OUTCAR_AMIX=0.1_AMIX_MAG=0.4.

I would greatly appreciate any suggestions or insights on how to address this issue with convergence, particularly concerning the sensitivity of the magnetic ground state and the convergence parameters.

Thank you in advance for your assistance!

Dear Neeta Bisht,

this is indeed strange and unexpected behavior. Some initial ideas:

• Is it possible that you used the free energy TOTEN to compare, instead of checking the extrapolation for sigma -> 0? To avoid dependencies on SIGMA, it is recommended to run static calculations with ISMEAR = -5 for accurate total energies, where SIGMA is ignored.

• Please check if the final magnetic moments differ when you vary MAGMOM or the mixing parameters. As you indicated, the energy differences you see point to different minima in the potential energy landscape.

• Make sure that you are running your comparison calculations statically without relaxing the substrate to be sure that the differences in energy are not coming from slightly different positions. In your substrate INCAR, EDIFFG = -0.04, which might be a non-negligable force in your case.

After looking at your OUTCAR for the substrate alone, I noticed that the forces in the bottom fixed two layers are very high. The bottom layer wants to relax outwards, and the middle one in the opposite direction. Are you sure that the interlayer distance is reasonable (bulk-like)? The layer distance between bottom and middle is 0.12 Angstrom lower than between top and middle. This might be fine if you want to force a bulk-like layer distance for the bottom, but the large forces (1.3 eV/Angstrom for the Co atoms) trouble me a bit.

I will investigate a bit more on the substrate and see if I can reproduce your large energy differences.
Cheers, Michael

Dear Michael,
Thanks for your fast response.
Here are the answers to some of your queries on my calculation:
1. Is it possible that you used the free energy TOTEN to compare, instead of checking the extrapolation for sigma -> 0?
All the energies I mentioned in my previous post are actually the extrapolation for sigma -> 0. Thanks for the ISMEARE=-5 suggestion. As you suggested, I ran two follow-up single point calculation with ISMEAR=-5, one for my test calculation SIGMA = 0.02 eV: -627.97482516 eV and the other for the AMIX and AMIX_MAG refined parameters (TOTEN= -646.03206512 eV). The final scf energies did not change with ISMEAR=-5 for AMIX and AMIX_MAG refined parameters Job. But, for the SIGMA = 0.02 eV: -627.97482516 eV Job it resulst in the following error ( with TOTEN= -624.37265272 eV). :
|
| W W AA RRRRR N N II N N GGGG !!! |
| W W A A R R NN N II NN N G G !!! |
| W W A A R R N N N II N N N G !!! |
| W WW W AAAAAA RRRRR N N N II N N N G GGG ! |
| WW WW A A R R N NN II N NN G G |
| W W A A R R N N II N N GGGG !!! |
| |
| Tetrahedron method does not include variations of the Fermi |
| occupations, so forces and stress will be inaccurate. We suggest |
| using a different smearing scheme, ISMEAR = 1 or 2 for metals in |
| relaxations.

2. The final magnetic momet differs slightly in case of change in SIGMA but vary significantly in case when I use MAGMOM = 27*5.0 27*-5.0 54*0.0 instead of MAGMOM = 27*4.0 27*-4.0 54*0.0. I have attached the inputs files for these two jobs with folder name MAGMOM=4 and MAGMOM=5 Jobs. As you can see there that INCAR, POSCAR and POTCAR for these two jobs are identical.
3. The comparison calculations are indeed with the same POSCAR in that case one can also take CONTCAR because these are scf calculations and the positions are not changed that much. But as mentioned in my previous post, the comparison energy for different SIGMA doesn't vary that much. I found out that SIGMA doesn't affect the scf energy that much. But, I would like to mention that the INCAR for the substrate that I provided has EDIFFG=-0.04, however the comparison calculations are done for EDIFFG=-0.001 eV/A. Since, EDIFFG=-0.04 is not a strict criterion. I have attached an example INCAR (INCAR_comparison_caluclations) for my substrate job that yields -641.07897822 eV energy. I agree with you that the forces are very high on the bottom two layers, unrelaxed atoms. But, I was under the assumption that since they are F F F in the POSCAR, this basically means that the forces are calculated for those atoms, but as they are fixed and not relaxed, one could basically ignore those high forces because the bottom two layers are simply not relaxed. Please correct me if I am wrong, and do we really need to look at the forces on those atoms ? I have to check my optimized surface unit cell to be sure if I actually took care of the interlayer bulk distance and also some literature. But, I will get back to you on that. Thanks for running a reproducibility check on substrate energy. However, I would recommend using this INCAR_comparison_caluclations for running substrate calculations.

Regards
Neeta Bisht

Deer Neeta Bisht,

first of all, the warning you received is not an error. It is a warning not to use ISMEAR=-5 for metals. You have a bandgap in your system AND you are running no relaxation, so the warning can be safely ignored.

You are right that the forces on atoms fixed by selective dynamics are ignored in the optimization. However, if there are large forces, the positions are not optimal, and convergence might be worse. However, this is not relevant to the problem you are having with energy differences for different optimization settings.

I now ran a couple of static calculations myself, all of them with ISMEAR=-5. For all of them, I used the structure defined in your substrate/CONTCAR from your initial post:

1. MAGMOM = 4, ALGO = Fast

2. MAGMOM = 5, ALGO = Fast

3. MAGMOM = 4, ALGO = Fast, AMIX = 0.2, BMIX=0.0001, AMIX_MAG=0.8, BMIX_MAG=0.0001

4. MAGMOM = 4, ALGO = Normal, AMIX = 0.2, BMIX=0.0001, AMIX_MAG=0.8, BMIX_MAG=0.0001

The rest of the INCAR is defined closely to what you used, and KPOINTS, and POTCAR are also equivalent to yours:

System = CoO-AFMII  

PREC = Accurate 
LASPH = .TRUE.
LREAL =Auto
NELM=300
ISMEAR = -5
ISPIN = 2
ENCUT = 450
LMAXMIX = 4
LORBIT = 11
EDIFF = 1E-6
LDAU = .TRUE.
LDAUTYPE =2
LDAUL = 2 -1 
LDAUU = 3 0  

LDAUJ = 0 0  

LDAUPRINT = 0
IVDW = 11
KPAR = 2
NCORE = 5

Calculations 1 and 2 gave me the same magnetic moments, but calculation 2 needed 72 instead of 53 electronic steps to converge. Energy without entropy is
equivalent as well: -642.0280 eV.

However, after switching to linear mixing for calculation 3, I observed changes in the oxygen magnetic moments, which are significantly larger for most atoms (\(\sim\pm0.1\mu_B\) instead of \(\sim\pm0.01\mu_B\)). The energy is roughly 1.5 eV higher at -640.6582 eV. Convergence is quicker (42 steps), but apparently, the code gets stuck in a different minimum with an energy difference of about 14meV/atom.

In the last calculation, the more stable blocked-Davidson algorithm for electronic minimization again finds the deeper minimum with the smaller O moments and the lower energy.

Overall, I think you are right. The system has a couple of shallow minima and you have to be careful with the settings, which is true for many magnetic systems.

To get accurate adsorption energies, I would perform 3 relaxations (slab, molecule, and combination) using settings that give you good convergence: Not too large initial moments, ISMEAR = 0 with a reasonable sized SIGMA, probably standard mixing with appropriate MAXMIX.

Following perform static calculations with ISMEAR = -5 and ALGO = Normal and use the energy values from those to compute the adsorption energy. Note that dipole corrections might play a role for the molecule and the combined system!

The OUTCARs of the calculations described above are attached, please do not hesitate to ask further questions!
Michael

Dear Michael,
Thanks for your calculations and your support.
As you mentioned:
To get accurate adsorption energies, I would perform 3 relaxations (slab, molecule, and combination) using settings that give you good convergence: Not too large initial moments, ISMEAR = 0 with a reasonable sized SIGMA, probably standard mixing with appropriate MAXMIX.

When you say not too large initial moments, does it mean the initial MAGMOM in the INCAR ? So, far I am dealing with only 4 and 5 as the initial MAGMOM value.

I looked into the OUTCAR_1, OUTCAR_2, OUTCAR_3, OUTCAR_4. I tried to reproduce your tests. Attached my OUTCAR_1, OUTCAR_2 in the folder called Reproduce_Michael_Test/without-KPAR. As you can see the results are quite different from yours, however, when I use NCORE=8 in combination with KPAR=2 in my INCAR (Reproduce_Michael_Test/KPAR), I am able to get the same energies as yours OUTCAR_1, OUTCAR_3. In either cases, the MAGMOM-5 and MAGMOM-4 does not come out to be same for me. Also, my OUTCAR_4 is scf unconverged in 300 iterations. I am rerunning this one with NELM=900. I have a query about your OUTCAR_1 and OUTCAR_2 results. I comapred the two OUTCAR's and found out one flag to be different i.e NELMDL. for OUTCAR_1 it is 0 and for OUTCAR_2 it is -5. NELMDL=0 indicates that WAVECAR is present. But, I think you started these calucltaions scartch ? Without reading the WAVECAR from the MAGMOM = 5, ALGO = Fast calulcations ? I am confused here, why these two flags are different. This could be the reason for the same energies for MAGMOM-4 and MAGMOM-5 in your case ? I ran the same scf run with two flags LDIPOL = .TRUE., IDIPOL = 3, the energies I got the sigma energy -641.35134931 eV.

Apart from this, I would like to mention that for all my comparison tests posted in my initial POST in the forum does not use substrate/CONTCAR because substrate CONTCAR was obatined with the following parameters (as given substarte/INCAR):

ALGO = Fast
PREC = Accurate
ADDGRID = .TRUE.
LASPH = .TRUE.
LREAL =Auto
AMIX=0.2
BMIX=0.0001
AMIX_MAG=0.8
BMIX_MAG=0.0001
NELM=300
ISMEAR = 0
SIGMA = 0.04
ISPIN = 2
ENCUT = 450

LMAXMIX = 4
MAGMOM = 27*5.0 27*-5.0 54*0.0

NSW =300
IBRION = 1
POTIM = 0.4
EDIFF = 1E-6
EDIFFG = -0.04
ISIF = 0

Since, here the EDIFFG criteion is not strict enough and the system is not that big, so I reoptimize the structure with substrate/POSCAR (substarte/POSCAR is formed from a optimized surface unit cell of CoO(100)) as my starting structure and with the following INCAR:
System = CoO-AFMII

ALGO = Fast
PREC = Accurate
ADDGRID = .TRUE.
LASPH = .TRUE.
LREAL =Auto
AMIX=0.2
BMIX=0.0001
AMIX_MAG=0.8
BMIX_MAG=0.0001
NELM=300
ISMEAR = 0
SIGMA = 0.01
ISPIN = 2
ENCUT = 450

LMAXMIX = 4
MAGMOM = 27*4.0 27*-4.0 54*0.0

NSW =300
IBRION = 1
POTIM = 0.4
EDIFF = 1E-6
EDIFFG = -0.001
ISIF = 0

LDAU = .TRUE.
LDAUTYPE =2
LDAUL = 2 -1
LDAUU = 3 0
LDAUJ = 0 0
LDAUPRINT = 0
IVDW = 11
NCORE = 8

The changes are essentially:
EDIFFG = -0.001 eV/A
SIGMA = 0.01 eV

This resulted in the energy -641.07897564 eV CONTCAR (attached file CONTCAR inside the folder refined_INCAR_optimization). All my mentioned comparsion single point calculations are done with this CONTACR. Additionally, if you run a scf calculations with the INCAR and CONTACR present in refined_INCAR_optimization which uses a reduce mixing parameter (AMIX=0.1 and AMIX_MAG=0.4 with BMIX=0.0001 and BMIX_MAG=0.0001) you will get a TOTEN = -646.03206512 eV. Which is even lower than all my previous tests. This energy value is also stable with ISMEAR=-5. I have done different test with combination of mixing parameter i.e. default as well as different combination of AMIX and AMIX_MAG also BMIX and BMIX_MAG. Every time I get something higher than -646.03206512 eV. The -646.03206512 eV run is present in the attached file folder JOB_AMIX=0.1_AMIX_MAG=0.4. It would be very helpful if you can look at this results and try to run the similar caluclations and check if this is reproducable. What is intresting here is that I can reproduce this lowest energy job results for SIGMA-0.01 and MAGMOM=5 as starting magnetic moment. So, far my guess is that the global minimum could be somewhat near to this value. Why I think that ? Because, if you run the ionic optimization with this Job INCAR and substarte/POSCAR as the starting structure you will get energy: -646.56425466 eV.

I appreciate your help in this matter.

Regards
Neeta Bisht

Dear Neeta,

When you say not too large initial moments, does it mean the initial MAGMOM in the INCAR ? So, far I am dealing with only 4 and 5 as the initial MAGMOM value.

Exactly. The MAGMOM tag sets the initial magnetic moments, and they should be selected higher than the expected moments, but not too high. We recommend 1.2 to 1.5 times the experimental value. We know that we should arrive at around 2.7 for the Co atoms in this system, so 2.7*1.5=4.05 is the upper limit I would try. And indeed, convergence is better for 4.0 than 5.0

I looked into the OUTCAR_1, OUTCAR_2, OUTCAR_3, OUTCAR_4. I tried to reproduce your tests. Attached my OUTCAR_1, OUTCAR_2 in the folder called Reproduce_Michael_Test/without-KPAR. As you can see the results are quite different from yours, however, when I use NCORE=8 in combination with KPAR=2 in my INCAR (Reproduce_Michael_Test/KPAR), I am able to get the same energies as yours OUTCAR_1, OUTCAR_3. In either cases, the MAGMOM-5 and MAGMOM-4 does not come out to be same for me. Also, my OUTCAR_4 is scf unconverged in 300 iterations. I am rerunning this one with NELM=900.

Apparently convergence is quite depended on the number of bands, NBANDS, as well for this system. I have 486 bands in all my calculations, you have 520 for KPAR = 2 and 494 for KPAR = 0. Note that your NCORE parameter is reset internally to 1, since your number of mpi-ranks/KPAR = 52, which is not devisible by 8. This is most certainly not ideal. Consult the wiki pages about parallelization, KPAR, and NCORE to understand better how to set these flags.
If your calculation does not converge within 300 steps, I would rather change some settings, than increase the number of steps.

I comapred the two OUTCAR's and found out one flag to be different i.e NELMDL. for OUTCAR_1 it is 0 and for OUTCAR_2 it is -5. NELMDL=0 indicates that WAVECAR is present. But, I think you started these calucltaions scartch ? Without reading the WAVECAR from the MAGMOM = 5, ALGO = Fast calulcations ? I am confused here, why these two flags are different.

I am very sorry for this confusion. I should have commented on this in my previous post. Indeed both calculations are done from scratch, but I restarted the first one after killing it (because I adjusted the parallelization settings I think). I did not delete the WAVECAR, which, just after the calculation started was empty. However, since the file was present, NELMDL is 0 in the OUTCAR. If you search for ISTART however, which indicates if the WAVECAR is used to initialize the orbitals, you will see that it is 0 for both calculations. And NELMDL is reset to the default (-5) if the WAVECAR is not read. So the calculations are comparable.

Additionally, if you run a scf calculations with the INCAR and CONTACR present in refined_INCAR_optimization which uses a reduce mixing parameter (AMIX=0.1 and AMIX_MAG=0.4 with BMIX=0.0001 and BMIX_MAG=0.0001) you will get a TOTEN = -646.03206512 eV. Which is even lower than all my previous tests. This energy value is also stable with ISMEAR=-5. I have done different test with combination of mixing parameter i.e. default as well as different combination of AMIX and AMIX_MAG also BMIX and BMIX_MAG. Every time I get something higher than -646.03206512 eV. The -646.03206512 eV run is present in the attached file folder JOB_AMIX=0.1_AMIX_MAG=0.4. It would be very helpful if you can look at this results and try to run the similar caluclations and check if this is reproducable. What is intresting here is that I can reproduce this lowest energy job results for SIGMA-0.01 and MAGMOM=5 as starting magnetic moment. So, far my guess is that the global minimum could be somewhat near to this value. Why I think that ? Because, if you run the ionic optimization with this Job INCAR and substarte/POSCAR as the starting structure you will get energy: -646.56425466 eV.

It sounds like you have found your minimum then. The only thing that seems slightly off to me is that some of your O atoms have a higher moment than I would expect for a low-energy solution. I will run some calculations today and try to confirm.

I still want you to double-check the converged bulk-lattice constant of rocksalt CoO with your employed AFM structure and computational settings (POTCAR choice, plane wave energy cutoff, and k-spacing) to see if your bottom two layers are indeed correctly spaced.
Note that you will likely make larger errors in the adsorption energy due to slab-thickness and structure, vdW approximation, and choice of density functional than you will by searching for the last 50meV (<0.5 meV/atom!) in the magnetic system.

I will get back to you as soon as my calculations are done,
Michael

Dear Michael,
Thank you very much on your constant feedback on my calculation. I will definitely read about the careful selection of NCORE, KPAR and its effect on NBANDS. I absolutely agree with you that the bulk CoO needs to be optimized properly. Hence, I started reading papers with DFT and experimental studies on CoO. I found three important references here:
1. https://journals.aps.org/prb/pdf/10.110 ... .75.104306 - Bulk lattice parameter 4.26 A
2. https://www.wellesu.com/10.1063/1.4937009 - Bulk lattice parameter 4.27 A
3. https://www.wellesu.com/10.1080/14328917.2017.1414110 - 4.26 A
4. DOI: https://doi.org/10.30880/ijie.2019.11.07.030 - 4.27 A

Ref1 is good but they worked with Rhombohedral (slightly distorted along 111 plane compare to CoO cubic ) for AFM-II. I have observed that for simplicity and computational efficiency, people often work with cubic structure (although KS-DFT is 0 K method and in that case CoO below Neel temperature have a slightly distorted - Rhombohedral phase). Nevertheless, I think, the cubic face with imposed AFM II ordering in INCAR is a good enough for adsorption energy calculations.
In the current calculations, I worked with a bulk unit cell with 4.25 Angstrom. I am now trying to look into my bulk optimisation of CoO bulk and I will get back to you. This is a crucial factor indeed. Do you think it is worth the time to look at the Energy convergence graphs and PES for Bulk CoO ? This will be fast.

Regards
Neeta Bisht

Dear Neeta Bisht,

I encourage you to closely study the bulk system before you explore surfaces and thin films. If it is OK to ignore the rhombohedral distortion is a scientific question that you have to decide on your own or within your research group.

Regarding my calculations on the relaxed slab:
I ran two, one with the standard mixing parameters, and one with linear mixing, which gives you the lowest energy (-646.032eV)
On my machine, running on 20 cores with KPAR = 2 and NCORE = 5 with 486 bands, the situation is quite different: I get -641.352 eV for the standard mixing, and -627.378eV with the linear mixing, which also converges much slower. Again, the magnetic subsystem is responsible. The average absolute Co moment for standard mixing is 2.65, while it is only 2.53 for linear mixing. The average absolute O moment is also higher for the standard mixing, with 0.09 instead of 0.06.

Note that my previous lower energy, for the unrelaxed slab, is actually lower than for the relaxed slab, at -642.038eV. The mean absolute Co moment is slightly higher at 2.70, while the mean absolute O moment stays small at 0.06. This means that the relaxed positions might also depend on the number of bands and how the magnetic subsystem evolves during the relaxation.

From all of your efforts and also my own couple of calculations, I conclude that indeed the PES of the slab is very flat with quite a few local minima. I suggest you monitor the magnetic moments in all your calculations carefully when computing the adsorption energy.

Please let me know if your initial slab setup was OK with respect to your bulk studies. Maybe adding an additional layer of Co100 would also help stabilize the magnetic system, but this is hard to predict.

Cheers, Michael

Dear Michael,
Thanks for your calculations and support again. Over the weekend, I looked into the bulk system carefully.
I have carried out a systematic bulk optimisation of CoO to get a suitable value for the lattice parameter and bulk inter-layer distance.
In my previous calculation, I used the AFM-II CoO(100) slab structure built from the non optimised bulk CoO. The lattice parameter of the non optimised CoO100 slab was 4.25 Angstrom with an inter-layer distance of 2.125 Angstrom. The experimental lattice parameter is 4.26-4.27 Angstrom (https://journals.aps.org/prb/pdf/10.110 ... .75.104306).

However, I obtained a 2D PES for my CoO(100) surface with 3 layers using a inter-layer distance between bottom two layer 1.98 Angstrom which is 0.15 Angstrom less than the experimental bulk inter-layer distance (Exp- 2.13 Angstrom).

I tried to converge the lattice parameter of the bulk CoO using GGA-PBE. Please have a look at the attached plots:

Plot_1.png - Lattice constant a, dependence upon U for GGA+U calculations. Its shows that U=2.6 eV and J=1 reproduce the experimentally calculated 4.26 Angstrom lattice parameter.
Plot_2.png - PES for U=3.0 eV J=0 eV, MAGMOM=2*5.0 2*-5.0 4*0.0 - Minima at 4.30 Angstrom (0.04 Angstrom Deviation)
Plot_3.png - MAGMOM dependence in U parameter. Here, J=1 in all the calculations. The shaded red region is the reported theoretical values.

When I run a calculations with U=2.6 eV and J=0 which mean taking U value 2.6 eV instead of 3.0 eV in Plot_2.png calculations . I get the same lattice minima as 1. i.e. 4.30 Angstrom.
It seems like J=1 in combination with U=2.5 and 2.6 gives me the experimentally obtained lattice constant value. In that case while running again the substrate with the fixed bulk inter-layer distance for bottom two layers, should I consider the same U and J parameter for the substrate ?

Regards
Neeta Bisht

Dear Neeta Bisht,

This is getting a bit off-topic I fear. Which values of U and J (mind that only U-J is relevant in Dudarevs approach that you selected with LDAUTYPE = 2) is a different problem and depends on what you want to focus on. Generally, it will not be possible to optimize for the size of the magnetic moments and for that lattice parameter at the same time. It is a tradeoff you have to think about and communicate when you present your results.

You should also be careful when optimizing the distance between the fixed layers and the lateral lattice parameter. I assume you fixed them because you want to represent the bulk on the bottom of your slab and only have the top relax to model a real surface on that side. So having strain on the bottom two layers is a good thing probably. But the forces in your initial calculations were large, so I was double-checking that you thought about what lattice parameter to use. The experimental lattice parameter, or the relaxed bulk lattice parameter for your selected U-J, are both reasonable choices to make. I suggest to discuss this with your supervisor or co-workers.

I also want to state again that you will have to add dipole corrections when you adsorb the molecule! If the top (relaxed) layer of your slab is buckling a lot, or if the interlayer distances are very different between the top and the middle layer and the middle and the bottom layer, they might be relevant for the bare slab as well.

If you still see a dependence on the smearing, you can use a 2x2x1 k-point mesh and ISMEAR = -5 again as well.

I am very happy to help with convergence issues, or with how to set VASP-specific parameters, but choosing U values or lattice parameters are scientific questions that are a bit out of the scope of this forum.

But I do not want to discourage you from asking more questions on VASP usage in this thread at all and I am happy to help more,
Cheers, Michael