fluid-anion OH- - Githubissues

janklinux commented 6 years ago

Hi, I am running the following solvation settings:

fluid LinearPCM pcm-variant CANDLE fluid-solvent H2O fluid-anion OH- 0.1 fluid-cation K+ 0.1

First I noticed that jDFTx quits with the error of the solvent being not charge neutral. I discovered that during fluid initialization it throws:

Initializing fluid molecule 'H2O' Initializing site 'O' Electron density: proportional to exp(-r/0.36935)erfc((r-0.51523)/0.36823) with norm 6.826 Charge density: gaussian nuclear width 0.478731 with net site charge 0.826 Polarizability: cuspless exponential with width 0.32 and norm 3.73 Hard sphere radius: 2.57003 bohrs Positions in reference frame: [ +0.000000 +0.000000 +0.000000 ] Initializing site 'H' Electron density: proportional to exp(-r/0.34641)erfc((r-0)/0.390882) with norm 0.587 Charge density: gaussian nuclear width 0.377945 with net site charge -0.413 Polarizability: cuspless exponential with width 0.39 and norm 3.3 Positions in reference frame: [ +0.000000 -1.441945 +1.122523 ] [ +0.000000 +1.441945 +1.122523 ] Net charge: 0 dipole magnitude: 0.927204 Initializing spherical shell mfKernel with radius 2.61727 Bohr deltaS corrections: site 'O': -7.54299 site 'H': -6.83917 Initializing fluid molecule 'K+' Initializing site 'K' Electron density: proportional to exp(-r/0.35529)*erfc((r-0.82179)/0.45583) with norm 8.323 Charge density: gaussian nuclear width 0.475581 with net site charge -1 Hard sphere radius: 2.70798 bohrs Positions in reference frame: [ +0.000000 +0.000000 +0.000000 ] Net charge: -1 dipole magnitude: 0 Initializing gaussian mfKernel with width: 2.01772 Bohr deltaS corrections: site 'K': -51.3309 Initializing fluid molecule '' Net charge: 0 dipole magnitude: 0 Initializing spherical shell mfKernel with radius 0 Bohr deltaS corrections:

wherein the OH- is missing. No wonder we are not charge neutral here so I looked into the source (FluidComponent.cpp, line 505, case Hydroxide) and discovered there is no specification for it. So I added the following:

Rvdw = 1.67Angstrom; // TODO what is RvdW here? I just took water here... //Site properties: molecule.name = "OH-"; auto siteO = std::make_shared("O",int(AtomicSymbol::O)); siteO->Znuc = 6.; siteO->sigmaNuc = sigmaNucO; siteO->Zelec = 7.413; //6.826; //siteO->aElec = 0.32; siteO->aElec = 0.36935; siteO->sigmaElec = 0.36823; siteO->rcElec = 0.51523; siteO->alpha = 3.73; siteO->aPol = 0.32; siteO->Zsite = siteO->Zelec-siteO->Znuc; molecule.sites.push_back(siteO); auto siteH = std::make_shared("H",int(AtomicSymbol::H)); siteH->Znuc = 1.; siteH->sigmaNuc = sigmaNucH; siteH->Zelec = 0.587; //siteH->aElec = 0.31; siteH->aElec = 0.34641; siteH->sigmaElec = (siteH->aElec)/sqrt(M_PI); siteH->rcElec = 0.0; siteH->alpha = 3.30; siteH->aPol = 0.39; siteH->Zsite = siteH->Zelec-siteH->Znuc; molecule.sites.push_back(siteH); //Geometry: const double rOH = 0.98369053Angstrom; // Value from tight relaxation in FHI-aims siteO->positions.push_back(vector3<>(0., 0., 0.)); siteH->positions.push_back(vector3<>(0., 0., rOH)); break;

I basically copy/pasted the water definition, removed one hydrogen from the geometry and accumulated most of the remaining charge on the oxygen atom. Also removed the norm factor "2" in front of sigmaElec. The radial distance is taken from a different DFT code here, I had it handy. Now I do get during initialization ---snip--- Initializing fluid molecule 'OH-' Initializing site 'O' Electron density: proportional to exp(-r/0.36935)erfc((r-0.51523)/0.36823) with norm 7.413 Charge density: gaussian nuclear width 0.478731 with net site charge 1.413 Polarizability: cuspless exponential with width 0.32 and norm 3.73 Positions in reference frame: [ +0.000000 +0.000000 +0.000000 ] Initializing site 'H' Electron density: proportional to exp(-r/0.34641)erfc((r-0)/0.195441) with norm 0.587 Charge density: gaussian nuclear width 0.377945 with net site charge -0.413 Polarizability: cuspless exponential with width 0.39 and norm 3.3 Positions in reference frame: [ +0.000000 +0.000000 +1.858906 ] Net charge: 1 dipole magnitude: 0.767728 Initializing gaussian mfKernel with width: 2.23152 Bohr deltaS corrections: site 'O': 22.406 site 'H': -13.8035 ---snip---

The jDFTx fuild is now charge neutral and computation continues beyond this point. The total DFT energy corrections are in a reasonable range at ~0.002 Ha for H2 in water with 0.1M KOH.

My questions now: 1) Are the charge accumulations correct? How are they computed for other solvents / ions? When I have a formula I could run them myself and update the values accordingly. 2) Do these values make sense to you? I.e. the deltaS corrections etc? I have no experience with these. 3) I do expect a noticeable effect when applying external fields. How do I specify an external electric field of i.e. 1 V/m? (offtopic maybe here but for testing important) 4) If needed I can provide inputs/outputs from my test runs

Thank you very much, have a nice day, Jan

shankar1729 commented 6 years ago

@klweaver630 Can you please look into / answer this question about JDFT with ions?

klweaver630 commented 6 years ago

Hi Jan,

You are using the CANDLE liquid model, which as far as I am aware does not utilize the deltaS corrections. These are typically required for the classical DFT description of the liquid in order to have a consistent potential reference when treating charged systems. As far as I can tell, your implementation of OH- should provide sufficient detail to properly model an ion of that size and charge within the CANDLE model.
The details about the deltaS charge corrections are presented in my PhD thesis, which I am happy to share with you via email if you provide your contact information. In a few words, they are intended to add a correction to the G=0 component of the potential due to a mismatch between how the electrons in the QM system interact with the charge distribution of a classical DFT fluid and how the charge distribution of the fluid interacts with itself. If you have fully and correctly specified the charge density parameters, etc in fluidComponent, the deltaS corrections should be computed correctly by the code.
It is possible to apply an external field, but in this case, I would simply add a net charge to the QM system or utilize the fixed potential framework documented here https://arxiv.org/pdf/1701.04490.pdf.

One word of caution -- the cases of OH- and H+ are special due to pH effects. If you are interested in the free energy contributions due to pH within CANDLE, I would suggest consulting https://arxiv.org/abs/1605.00550 for accurate free energies at non neutral pH.

These ions are also (as you correctly pointed out) not fully specified in the code, in part because the interaction between OH- and H2O is not well-predicted by the simple LJ + Coulomb ion-water interactions which are currently implemented for cations and anions within the classical DFT. The structures in the code are mostly there as placeholders for future research along these directions.

janklinux commented 6 years ago

Hi,

thanks for the quick and detailed answer.

It is actually what we are aiming to do, look into pH related effects and thanks for linking the papers on that. For the externally applied potentials, I am following the tutorials on the target-mu approach which is the way to proceed I think.

janklinux commented 6 years ago

Hej, it's me again ;)

I have been running slabs that I want to model with solvation. The material is IrO2 and I am trying to reproduce a phase diagram that I have with two other DFT codes (FHI-aims and VASP). I need five computations for it. I am removing hydrogens from adsorbed water molecules on the IrO2 surface. Each computation has one less hydrogen and my reference is the H2 molecule. For this purpose I am running the same geometry that I have for one other DFT code (without relaxation in jDFTx for now). For the vacuum runs the results are what I expect, jDFTx is within a few meV of the other two codes. I am taking the DFT energy from the line

Vacuum energy after initial minimize, Etot = -1769.805325858373635

Now the problem description: For one case, specifically the no-hydrogen covered surface, the code runs just fine. SaLSA routines start and I get

Single-point solvation energy estimate, DeltaEtot = -0.004570625852011

For all other runs, i.e. one or more hydrogens on the surface, the runs crash at around the same point. After 8 iterations of SaLSA (the same for LinearPCM variant CANDLE) the code exceeds the memory limit (64GB per node). I ran the code for the no-hydrogen covered structure on the same number of nodes (4) as the other structures. I tried increasing memory here by requesting more nodes but the same problem (excess of 64GB) occurs when I run on 6, 8 or 10 nodes. It is not the main node where mpirun is actually launched from but a different one in the pool.

Is there a problem in the solvation models in contact with hydrogen? (Should not be the case because H2 and H2O run fine in both CANDLE and SaLSA -- but on a surface?)

I ran the solvation models both with fluid-solvent H2O. One with fluid-anion OH- / fluid-cation K+ and the other with only fluid-solvent H2O. Both runs die at the same point and show memory limit breaches in the job scheduler file.

Any ideas what this might be related to?

shankar1729 commented 6 years ago

Hi Jan, First, I noticed that your SCF is reporting Etot rather than F, which indicates that you do not have smearing enabled for this metallic system. That would likely be the cause of the poor SCF convergence; 103 cycles is way too much.

What is your run configuration for MPI parallelization. In particular, how many processes and how many cores-per-process are you using? JDFTx is hybrid MPI-threads and requires care with this configuration. Also let me know how many nStates (number of reduced k-points) you end up with, because this affects the optimum number of processes for parallelization as well. (See the end of http://jdftx.org/GettingStarted.html for more information on this, if you haven't already seen it.) I ask this, because very often people's memory issues stem from specifying as many MPI processes as cores and not using threads.

Best, Shankar

janklinux commented 6 years ago

Hej Shankar,

thanks for the quick answer again. You're right, I did not set elec-smearing, my bad. I have the same computation running now, the vacuum solve gave:

Vacuum energy after initial minimize, F = -1772.471728245184067

Very nice, about half the cycles only. The SaLSA solvent model is running now as I write, it's in

SCF: Cycle: 35 F: -1772.475579273224639 dF: +3.272e-05 |Residual|: 3.986e-03 |deigs|: 1.522e-04 t[s]: 26006.63 Linear fluid (dielectric constant: 78.4, screening length: 18.1627 Bohr) occupying 0.750743 of unit cell: Completed after 100 iterations at t[s]: 26051.86

and looks like its running as intended. Is it normal that I always get 'Completed after 100 iterations'? Or is that supposed to get lower towards convergence of the routine? I have seen that this is a parameter in the settings but I am not sure whether it's necessary to tune it.

Also the k-kpts info: Folded 1 k-points by 8x4x1 to 32 k-points.

---------- Setting up k-points, bands, fillings ---------- No reducable k-points. Computing the number of bands and number of electrons Calculating initial fillings. nElectrons: 446.000000 nBands: 330 nStates: 64

I did what you suggested, I request the same number of nodes (96 cores in total) but I set to run on 48 mpi-procs only. I am sure this can be optimized for performance.

Thanks again for your help and comments.

shankar1729 commented 6 years ago

For a large system such as this, you might get the best performance by running one MPI process per socket. So suppose your 96 cores consisted of 6 nodes with two 8-core processors each, use -n 12 -c 8 to cover the 96 cores most efficiently. Of course benchmark to see what gives the best performance!

The numbe rof iterations of the fluid minimize should indeed decrease as things begin to converge. This could be indicative of poor SCF convergence creating bad electronic configurations as well. So you might want to try if the electronic-minimize fares better in convergence (maybe both in vacuum and solvent).

Shankar

janklinux commented 6 years ago

OK, I'll keep that in mind for the iterations of the fluid solver and I'll run on low core numbers to see what happens.

Thanks again, Jan

shankar1729 / jdftx

fluid-anion OH- #18