Forces with Soft-sphere model

[BCAyers2000]

Dear Environ Development Team,

I hope this message finds you well. I'd like to express my appreciation for the memory issue fix in the update from Environ 3.0 to 3.1. This improvement is significant for our work.

However, I've encountered a new challenge with force convergence when using soft-spheres in version 3.1. Previously, with version 3.0, we could consistently achieve convergence to 1E-4 Ry/Bohr. Now, using 3.1, the calculation stalls at a gradient of 1.5 Ry/Bohr after approximately 11 iterations and remains fixed at this point for over 65 iterations before I had to terminate the process. I've thoroughly checked our setup, recompiled the software, and validated this behavior on multiple HPC systems. As nothing else has changed in our workflow besides the Environ update, I suspect this issue may be related to the recent modifications.

For your reference, I've attached our environ.in file and the truncated pw.out file: files.gz

This convergence issue is critical for us, as we need to run simulations on larger systems. While the memory fix allows us to attempt these larger calculations, the lack of force convergence prevents us from obtaining meaningful results. Any insights or assistance you could provide would be greatly appreciated. Thank you for your time and continued support.

Best regards, Brad Ayers

[eric-read]

Can you send your QE input and environ.debug files?

[olivieroandreussi]

Thank you very much Brad for reporting this issue. We released 3.1 in a bit of a hurry due to an event we hosted last week. I am sorry we should have checked the forces of the soft-sphere version, indeed the clean up of the old code was not complete. I have committed a cleaner version that passes all our internal tests related to soft-spheres. Would you mind testing it and let us know if you have issues? We will finish testing the other parts and will release 3.1.1 if everything looks good. Oliviero

[BCAyers2000]

Thank you so much for the quick response!

I’ll get started on compiling now and let you know if everything works out. I’ve also attached the requested files below, just in case they’re helpful. They’re named Environ_3.0 and 3.1, and each contains the environ.in, pw.in, pw.out, and environ.debug files.

environ_debugging.tar.gz

Also I know this is not the best forum for this, but when using the gcscf implementation, certain structures fail to get good forces even with old environ, would you be able to look over my input file and advise as I'm not sure what else I can try, I have done solvent-aware, solvent distance , electrolyte distance etc.

Best, Brad

[eric-read]

I believe that discussion would be better suited on the Environ google group

[BCAyers2000]

Hi Both,

After your fix the forces seems to work fine now!

kind regards, Brad

[BCAyers2000]

Hi All ,

something new that I have noticed is the emergence of this error message at the start of every bfgs iteration :

Environ Warning (1005): wrong integral of erfc function

Environ Warning (1005): wrong integral of erfc function

Environ Warning (1005): wrong integral of erfc function

Environ Warning (1005): wrong integral of erfc function

It repeats this quite a few times as well , never during the scf but at the start it is common, I am not sure this happened often before.

[olivieroandreussi]

Thanks for the feedback. This usually happens when either the cell is small in at least one direction (i.e. if you are doing a periodic slab with a 1x1 unit cell in the xy plane) or if your grid is quite coarse (i.e. if ecutrho is very low, say in the 100 Ry range, instead of 300)

[olivieroandreussi]

Are you specifying any of the soft-sphere parameters by hand? say the spread? This error may also happen if you define not-so-soft soft spheres, even if you use fine grids.

[BCAyers2000]

Hiya,

I'm working with this input file for Environ:

&ENVIRON verbose = 1 cion(1) = 1 cion(2) = 1 zion(1) = 1 zion(2) = -1 rion = 3.D0
system_dim = 2 system_axis = 3 environ_thr = 1.d-1
env_pressure = 0.0 temperature = 298.15 environ_type = 'input' env_electrolyte_ntyp = 2 env_surface_tension = 37.3 env_static_permittivity = 90.7 / &BOUNDARY alpha = 1.0 radius_mode ='muff' solvent_mode = 'ionic' solvationrad(1) = 4.01 electrolyte_mode = 'ionic' electrolyte_alpha = 1.0D0 / &ELECTROSTATIC pbc_axis=3 pbc_dim =2 tol=1.0e-10 problem='modpb' inner_tol = 1.0e-15 pbc_correction = 'parabolic' /

The default Li radii in Environ weren't working for me - the SCF just wouldn't converge. So I dug around in some other softwares (ONETEP) and found they use values from Alvarez's "A cartography of the van der Waals territories". Using those seems to work much better.

My cell block looks like this:

CELL_PARAMETERS angstrom 4.85077600855226 0.00000000000000 0.25252131244487 2.40566940261534 8.04605936887865 0.50504262488974 0.00000000000000 0.00000000000000 83.86295788450576

With ecutrho at 320, I think this should just fit two soft-spheres without them interacting across periodic boundaries.

While I've got the chance , I'd like to as a few more questions if that's ok:

1. I'm trying to model Li+ and PF6- ions in an electrolyte using modpb. I've set rion to 3 Bohr, which is the average of Li+ (1.341 Bohr) and PF6- (4.8 Bohr). Does that make sense, or is there a better approach?

2. For the SSCS method with GC-DFT, would turning on field-aware be helpful?

Thanks for any help you can give me!

Cheers, Brad

[olivieroandreussi]

Thanks for sharing your experience with Li radii, it sounds reasonable and useful to others. I may add this value to our muff table (modified uff). The cell seems large enough to have two spheres and the cutoff should be ok, I am not sure which sphere is giving the error... The fact that your cell is not orthorhombic and you are using the parabolic correction seems a bit dangerous to me, I wonder if this has anything to do with the warning message. On the other questions, 1. I think anything in the right ballpark should be fine, although I would probably go higher (Li+ is probably solvated and has a larger effective radius). I would probably try a check a posteriori to see how much a different value affects the results, hopefully not too much. 2. field-aware requires extra parameters and at this time we don't have a final recipe for those (Eric, who answered some questions above, is also working on this). I would also expect field-aware to make convergence more challenging, it is meant to improve accuracy in solvation free energies, not stability of electrochemical simulations. I have not played with GC-DFT enough to know the right tricks to help convergence, I believe things that help converging the SCF may help with the GC algorithm (reduce mixing_beta, add smearing, etc.). But I have not tried these myself. If you find anything that helps, let us know, as we are also starting to test GC more.

[BCAyers2000]

Thank you for your fast and detailed reply!

regarding using the parabolic correction, what would you recommend as using 'none' with modpb seems not to be supported?

kind regards, Brad

[olivieroandreussi]

Is there a way you can orthogonalize the cell? The components of a1 and a2 that are along z are small, I am not sure if those make sense (are really needed) for a slab calculation. Is this because of the cell vectors of the starting bulk (monoclinic or triclinic phase)? Or did you relax the cell with the vacuum and got this? In the latter case, I would try to do a more fancy cell relax, using celldofree and keeping the symmetry of the cell fixed (e.g. if you want to keep it hexagonal in the xy plane).

In any case, you are right, to do modpb you need the parabolic correction. It is also possible that the correction works fine for this kind of cell, but I believe in some parts of the implementation I assume a3 is orthogonal to the other two axis (e.g. the area perpendicular to z is computed by taking the volume and dividing it by the length of a3 along z). While some of these can be changed to a more general formulation, I am afraid that the parabolic correction is only well defined for a cell with a third axis perpendicular to the other two. This suggests I should add a warning to capture situations where this is not the case, I haven’t seen them before.

Oliviero

[BCAyers2000]

Hi again,

I hope this message finds you well. I apologise for the repeated inquiry, but I believe you're best positioned to clarify this matter. I'm seeking clarification on how Quantum ESPRESSO (QE) reports total charge, specifically whether it uses a negative or positive electron convention.

I've noticed some discrepancies between QE and Environ outputs:

1. Environ debug file reports final:

   • Number of electrons: 114
   • Total electronic charge: 114.3435725
   • Total charge: -5.6564275

2. QE output shows:

   • Total charge: 5.65642600 e

Given that my system started with 120 electrons, it appears that QE is reporting the charge correctly, indicating a loss of about 5.66 electrons. In contrast, Environ seems to calculate the charge as: Q = current_electrons - starting_electrons, hence the negative charge?

I'm particularly concerned because I've obtained negative surface energies when using GC-SCF and Environ together. The potential contribution is reported as:

pot.stat. contrib. (-muN) = -2.14218166 Ry

This is negative, despite having a positive charge and a negative chemical potential (μ = -5.15267399 eV). Logically, -μN should be positive in this case (-*- = +) and rather μN should be negative?

So, to cut a long ramble short do you know the following:

1. The charge convention used in QE output?
2. The correct interpretation of the potential contribution (-μN)?
3. Any insights on why I might be seeing negative surface energies?

Your expertise and clarification would be immensely helpful in deciphering these results. Thank you for your continued support and contributions to the scientific community. Your work is greatly appreciated!

Best regards, Brad

[eric-read]

Environ calculates the total charge by integration. However, if you are using an ionic interface, if the cell is too small in an axis such that the soft sphere of an atom overlaps with the soft spheres of its own periodic replicas, it can cause those integrations to become incorrect due to overlap.

[BCAyers2000]

Well It seems correct but rather a sign swap of QE, i.e. -5 instead of +5 e ?

my cell is :

CELL_PARAMETERS angstrom 6.86932225600000 0.00000000000000 0.00000000000000 -0.00000000000000 6.86932225600000 0.00000000000000 0.00000000000000 0.00000000000000 55.45597507600000

and my soft-sphere radii is 4.01 Bohr , so I don't think id get said interactions?

I just ran a quick constant charge calculation of +4e , and it does show up as -4e approximately in the output file so I guess it is the reverse in the environ debug file.

[eric-read]

I think Environ reports the number of electrons

[eric-read]

your cell parameters should be fine if your lattice vectors are larger than the soft sphere radii, although you do have to take into account the alpha scaling parameter

[eric-read]

Although as I have previously stated these conversations are generally better suited for the Environ google group