e: Identical Parameters Yield Different Outcomes in Multi-layer vs few-layer Systems

[Brad Ayers]

Dear Environ Development Team,

I am writing to report an intriguing convergence issue I've encountered whilst using Environ. I have two nearly identical systems - one with fewer layers that converges successfully, and another with additional layers that fails to converge, despite using identical environmental parameters.

What makes this particularly puzzling is that both the Polarization accuracy and Electrolyte accuracy converge normally in both cases. 

For reference, while quadrupole moments do show some differences between the systems (with the larger system showing values beyond Fortran's display limits), the primary concern is understanding why adding layers would prevent convergence despite well-behaved polarization and electrolyte parameters.

I have attached two tar.gz files containing the complete set of relevant files for both systems:

1. converged.tar.gz - Contains espresso.in, environ.in, espresso.out, and environ.debug for the working system
2. unconverged.tar.gz - Contains the same set of files for the non-converging system, stopped early as i have run this a few times and i know it doesn't converge. 

I would greatly appreciate your insights into what might be causing this behaviour and any suggestions for achieving convergence in the larger system.

Thank you for your time and assistance.

Best regards,
Brad

[Brad Ayers]

I wanted to follow up on this issue—adjusting the interface spread to 0.60 Bohr seems to resolve the convergence problem observed in the larger system. With this modification, the transition remains well-behaved while still reflecting a reasonable solvation response.

The updated parameters are as follows:

• Bondi radius of Li: 1.82 Bohr
• Scaling factor: 1.32
• Effective radius: 2.40 Bohr
• Interface spread (softness): 0.60 Bohr

From my tests, this adjustment ensures that both systems converge successfully while maintaining physically meaningful behaviour.

[andr...@gmail.com]

Thanks Brad for the updates and for the interesting question. I have not seen this before, but I know Li can be a challenging system. The solution you posted, increasing the spread, seems to suggest that the problem is with how close the interfaces get to the electrons of the system. It seems reasonable that pushing the interfaces away improves the stability of the calculation. 

My main concern with adjusting the spread is that, in principles, you need to retune the other solvation parameters for the new setup. In the past a student of mine tested using a spread of 0.8 and showed that we had to refit alpha to 1.21 to get the same accuracy on solvation free energies of neutral molecules. I see you are using specific radii and scaling factors, so I am guessing you chose them to reproduce some property you are interested in. I would make sure that those properties still make sense with the new spread.

An alternative way to get the same effect (pushing the boundary away from the electrons) is to use the solvent-aware correction, specifying a solvent radius in the BOUNDARY namelist. I believe this forum has other conversations where I mention this with more details. The solvent-aware was shown to help the convergence of another tricky system (TiO2) possibly for the same reasons you see in your applications. The advantage of using the solvent-aware correction is that it acts on specific parts of the interface, while leaving the parts that are more exposed unaffected. This means that in most cases the effect is stronger on the tricky regions, but the overall solvation is still ok without a reparametrization. However, be aware that the solvent-aware correction may introduce some small errors in the forces that may not be as negligible as the standard ones from Environ (if this is the case, you may want to decrease the convergence threshold on the forces for the relax calculations). 

Thanks again,

Oliviero

Oliviero Andreussi (he/him)
--
Associate Professor
Department of Chemistry and Biochemistry
Boise State University
Office: SCNC 316
Email: oliviero...@boisestate.edu
Web: http://www.materialab.org

[Brad Ayers]

Thank you for your response. I'd like to clarify my current approach and challenges:

I am working on simulating lithium slabs in ethylene carbonate (EC) with 1M electrolyte at varying voltages. Due to difficulties implementing SCCS, I've been utilising SSCS with Bondi radii and a 1.32 scaling factor, following the Fissicaro methodology. Whilst I acknowledge these parameters were originally determined for neutral molecules, finding reliable values for lithium solvation has proven particularly challenging.

My primary objective is to analyse the evolution of surface energy in relation to solvation and voltage. Thus far, I've conducted this analysis using Bondi radii and 1.32 scaling for voltages . However, when investigating voltage effects on the adsorption energy of F onto Li—which required larger symmetric surfaces—I encountered the significant issues I discussed above. The observed spread (0.6)  actually alters the surface energy values unfavourably, suggesting greater stability in vacuum than in solvent, which contradicts expected behaviour... so it seems you were likely right haha. 

A key constraint with the solvent-aware system is the size requirement: ethylene carbonate's radius of approximately 5 bohrs necessitates a larger cell, which significantly increases computational costs... 

Most relevant literature I've found uses VASP with VASPSOl++, employing a cavity size parameter nc (https://pubs.acs.org/doi/10.1021/acs.jctc.5b00170); however, I'm uncertain how to translate this parameter to equivalent environ parameters.  I would greatly appreciate any guidance, as my research group has limited experience with solvent/Li work, so I am going into this fairly blind, so any help or ideas would be grand. 

many thanks,

Brad

[andr...@gmail.com]

I wrote a long reply on the browser, but it may have been lost instead of being posted... Sorry if this is a redundant message (and sorry if I am going to keep this one shorter).

The system is tricky, I can see why using a large radius would make the calculation too expensive. However, using no solvent aware corresponds to no radius, so you could try to use a small radius different from zero and it would still make some sense. The goal is mostly to improve convergence, it does not need to be a physical parameter, sort of like using a smearing to help converging the SCF for a non metallic system.

You could also try to see if solvent-aware (SA) helps with converging the SCCS calculations, I know they have problems with Li, but maybe SA would improve. If you try it, it would be interesting to know it. 

We cannot reproduce VASPSol++, the interfaces are built differently. However, we have multiple ways to control the Environ interfaces and it may be helpful to know which ingredient is causing the problem: dielectric, electrolyte, non-electrostatic (surface and volume). As of now dielectric and non-electrostatic are defined on the same interface, which makes it less than ideal to fix this kind of issues. However the electrolyte has its own interface and it may be possible to just smear or push away that interface, without affecting the main solvation. You could try to turn off the different ingredients one at the time and see if the calculations converge. 

Hope this helps, let me know if you need more details.

Best,

Oliviero

[Brad Ayers]

Thank you for your helpful response regarding convergence issues. I wanted to follow up with some specific results from my attempts to model lithium in ethylene carbonate:

1. I performed a systematic scan of solvation radii from 1.5-5.5 Bohr:

• Calculations failed to converge for radii < 2.9 Bohr
• Convergence begins at 2.9 Bohr with surface energy of -5.5764 J/m²
• Surface energy becomes positive around 3.7 Bohr (0.1710 J/m²)
• Converges to the VASPSOL++ result for Li at 4.3 bohr (0.462  J/m²)

2. I then tried implementing MUFF radii with alpha scaling:

• Even with alpha up to 1.40, the calculations still failed to converge

      3. I tried the solvent aware (SA) scheme from 1-4 bohrs and SCCS never converged... 

Given that convergence only starts at relatively large radii (2.9 Bohr), should I be adjusting the rhomax/rhomin parameters rather than just using the defaults?
Any guidance would be greatly appreciated, particularly on whether I should focus on optimising the electrostatic parameters (rhomax/rhomin) or if I should just use the 4.3 bohr parameter.


Sorry for bombarding you with questions,
Best,
Brad

[andr...@gmail.com]

Thanks for the update! While you could play with rhomax/rhomin, I would feel even less comfortable adjusting those. I wanted to clarify that MUFF stands for modified UFF, and the modification was in the nitrogen parameters and meant to improve the accuracy of solvation energies. If you know a specific element has the 'wrong' radius for solvation properties, there is nothing wrong to adjust just that radius. You can adjust the parameter in the fortran file that contains all the radii, or you can specify the solvationrad for the atom type that corresponds to Li. I cannot remember how the input parsing works, I am not sure if you can specify muff and then adjust a single radius, I feel this should be the way it works, but I may remember wrong (the alternative would be to specify all the radii in the input). Maybe this is already what you were doing, if this is the case, and you are only changing a specific radius, I think it would be reasonably easy to justify the approach. 

Best,

Oliviero