calculation of the grand-canonical electronic energy about vaspsol

[rhennig]

Hi Xinnan, The energy of the slab is used for the y-axis, which plots the grand-canonical free energy, F, relative to the value at zero applied potential. This set of calculations changes the electron count on the slab in the VASPsol DFT calculations. From the DFT calculation, we extract the difference between the Fermi level and the electrostatic potential inside the electrolyte. Multiplying the excess charge with this FERMI_SHIFT provides the extra energy term ΔE_ref that is subtracted from the energy of the slab. This is described in the caption and the text. Let us know if you have further questions. Best regards, Richard On Feb 4, 2022, at 3:30 PM, Xinnan Mao ***@***.******@***.***>> wrote: [External Email] Hello, I saw FIG. 2 in the VASPsol publication https://aip.scitation.org/doi/10.1063/1.5132354<https://urldefense.proofpoint.com/v2/url?u=https-3A__aip.scitation.org_doi_10.1063_1.5132354&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=XFhs398m2tUV4CPT-IO43vvQr0Cdxc7LpXYEfp-OoFtSz3-cDQsk3Uul0oCwackg&s=ietmCxYK_XrgvYbuA56MTXQ1Z_hUhoq3uODoHBjBbsw&e=>. Is there any mistake in my understanding: 1. The energy of the charged slab in VASP output is never used. 2. The correction add to the energy of the neutral slab is given by change of Q U relative to the neutral slab, where Q = net_charge U = E_vacuum_level(in LOCPOT) - EFermi(in OUTCAR) - FERMI_SHIFT(in stdout) — Reply to this email directly, view it on GitHub<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_henniggroup_VASPsol_issues_61&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=XFhs398m2tUV4CPT-IO43vvQr0Cdxc7LpXYEfp-OoFtSz3-cDQsk3Uul0oCwackg&s=v40_Uroaqnfdw4JVZDMZsubcEsO7NWH_ghIz_T2_Y0s&e=>, or unsubscribe<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_ACCOT2UU2UD33JY2XJFRXHTUZQZPLANCNFSM5NSXVHOA&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=XFhs398m2tUV4CPT-IO43vvQr0Cdxc7LpXYEfp-OoFtSz3-cDQsk3Uul0oCwackg&s=xTUODEgUOEV9XTJiMCdPeqqqr6F9-Ojq6pLyKeabOms&e=>. Triage notifications on the go with GitHub Mobile for iOS<https://urldefense.proofpoint.com/v2/url?u=https-3A__apps.apple.com_app_apple-2Dstore_id1477376905-3Fct-3Dnotification-2Demail-26mt-3D8-26pt-3D524675&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=XFhs398m2tUV4CPT-IO43vvQr0Cdxc7LpXYEfp-OoFtSz3-cDQsk3Uul0oCwackg&s=wjyyukblGZbZDF-ksLGAC_5fnt1-eB1KQK815_ht8eQ&e=> or Android<https://urldefense.proofpoint.com/v2/url?u=https-3A__play.google.com_store_apps_details-3Fid-3Dcom.github.android-26referrer-3Dutm-5Fcampaign-253Dnotification-2Demail-2526utm-5Fmedium-253Demail-2526utm-5Fsource-253Dgithub&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=XFhs398m2tUV4CPT-IO43vvQr0Cdxc7LpXYEfp-OoFtSz3-cDQsk3Uul0oCwackg&s=hg8pLNd30LcXTXEd9yq5pwKx8Dp9_-qnRCsMF-ZQahw&e=>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[xnmao]

Dear Prof. @rhennig ,

I have calculated a set of charged (2×2)-Pt(111) slab. The results are as follows:

	net_charge	potential_energy	FERMI_SHIFT
0	-0.20	-93.781	2.041
1	-0.16	-93.494	2.041
2	-0.12	-93.204	2.042
3	-0.08	-92.910	2.042
4	-0.04	-92.612	2.042
5	0.00	-92.312	2.041
6	0.04	-92.007	2.041
7	0.08	-91.700	2.042
8	0.12	-91.389	2.042
9	0.16	-91.075	2.041
10	0.20	-90.757	2.041

Is that grand-canonical free energy F = potential_energy - net_charge*FERMI_SHIFT? It seems a linear relationship between net_charge and F, no quadratic curve.

Here is the POSCAR:

Pt                                      
   1.00000000000000     
     5.6416157960000008    0.0000000000000000    0.0000000000000000
     2.8208078980000004    4.8857825977275668    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.9095400125127817
   Pt
    16
Selective dynamics
Direct
  0.0000000000000000  0.0000000000000000  0.3427243190829614   T   T   T
  0.5000000000000000  0.0000000000000000  0.3427243190829614   T   T   T
  0.0000000000000000  0.5000000000000000  0.3427243190829614   T   T   T
  0.5000000000000000  0.5000000000000000  0.3427243190829614   T   T   T
  0.1666666666666643  0.1666666666666643  0.4474388781586569   F   F   F
  0.6666666666666643  0.1666666666666643  0.4474388781586569   F   F   F
  0.1666666666666643  0.6666666666666643  0.4474388781586569   F   F   F
  0.6666666666666643  0.6666666666666643  0.4474388781586569   F   F   F
  0.8333333333333357  0.3333333333333357  0.5525611218413431   F   F   F
  0.3333333333333357  0.3333333333333357  0.5525611218413431   F   F   F
  0.8333333333333357  0.8333333333333357  0.5525611218413431   F   F   F
  0.3333333333333357  0.8333333333333357  0.5525611218413431   F   F   F
  0.0000000000000000  0.0000000000000000  0.6572756809170386   T   T   T
  0.5000000000000000  0.0000000000000000  0.6572756809170386   T   T   T
  0.0000000000000000  0.5000000000000000  0.6572756809170386   T   T   T
  0.5000000000000000  0.5000000000000000  0.6572756809170386   T   T   T

and INCAR

INCAR created by Atomic Simulation Environment
 ENCUT = 520.000000
 KSPACING = 0.300000
 NELECT = 160.080000 # net_charge = -0.08
 EB_K = 78.360000
 TAU = 0.000000
 LAMBDA_D_K = 3.040000
 ALGO = Fast
 GGA = PE
 LCHARG = .FALSE.
 LVHAR = .TRUE.
 LSOL = .TRUE.
 LREAL = Auto

Best regards,
Xinnan

[rhennig]

Dear Xinnan, You are missing the the Fermi level in your output. The electrostatic potential of the slab is given by the calculated Fermi level referenced by the Fermi Shift. Best regards, Richard — Richard G. Hennig Professor and Alumni Professor of Materials Science and Engineering Associate Director, Quantum Theory Project Department of Materials Science and Engineering, University of Florida Email: ***@***.******@***.***> Web: http://hennig.mse.ufl.edu On Feb 5, 2022, at 1:34 PM, Xinnan Mao ***@***.******@***.***>> wrote: [External Email] Dear Prof. @rhennig<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_rhennig&d=DwMFaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=CIoP38NMAXgBnnLsCELUobX_Pv7ewLJQLbNCmAz_GzOv1MUkwpdIa8DraR-2PBp4&s=_OwYWqrY66ymidDGBMX4cVEDskj_PRacQVokVgJjKQA&e=> , I have calculated a set of charged (2×2)-Pt(111) slab. The results are as follows: net_charge potential_energy FERMI_SHIFT 0 -0.20 -93.781 2.041 1 -0.16 -93.494 2.041 2 -0.12 -93.204 2.042 3 -0.08 -92.910 2.042 4 -0.04 -92.612 2.042 5 0.00 -92.312 2.041 6 0.04 -92.007 2.041 7 0.08 -91.700 2.042 8 0.12 -91.389 2.042 9 0.16 -91.075 2.041 10 0.20 -90.757 2.041 Is that grand-canonical free energy F = potential_energy - net_charge*FERMI_SHIFT? It seems a linear relationship between net_charge and F, no quadratic curve. Here is the POSCAR: Pt 1.00000000000000 5.6416157960000008 0.0000000000000000 0.0000000000000000 2.8208078980000004 4.8857825977275668 0.0000000000000000 0.0000000000000000 0.0000000000000000 21.9095400125127817 Pt 16 Selective dynamics Direct 0.0000000000000000 0.0000000000000000 0.3427243190829614 T T T 0.5000000000000000 0.0000000000000000 0.3427243190829614 T T T 0.0000000000000000 0.5000000000000000 0.3427243190829614 T T T 0.5000000000000000 0.5000000000000000 0.3427243190829614 T T T 0.1666666666666643 0.1666666666666643 0.4474388781586569 F F F 0.6666666666666643 0.1666666666666643 0.4474388781586569 F F F 0.1666666666666643 0.6666666666666643 0.4474388781586569 F F F 0.6666666666666643 0.6666666666666643 0.4474388781586569 F F F 0.8333333333333357 0.3333333333333357 0.5525611218413431 F F F 0.3333333333333357 0.3333333333333357 0.5525611218413431 F F F 0.8333333333333357 0.8333333333333357 0.5525611218413431 F F F 0.3333333333333357 0.8333333333333357 0.5525611218413431 F F F 0.0000000000000000 0.0000000000000000 0.6572756809170386 T T T 0.5000000000000000 0.0000000000000000 0.6572756809170386 T T T 0.0000000000000000 0.5000000000000000 0.6572756809170386 T T T 0.5000000000000000 0.5000000000000000 0.6572756809170386 T T T and INCAR INCAR created by Atomic Simulation Environment ENCUT = 520.000000 KSPACING = 0.300000 NELECT = 160.080000 # net_charge = -0.08 EB_K = 78.360000 TAU = 0.000000 LAMBDA_D_K = 3.040000 ALGO = Fast GGA = PE LCHARG = .FALSE. LVHAR = .TRUE. LSOL = .TRUE. LREAL = Auto Best regards, Xinnan — Reply to this email directly, view it on GitHub<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_henniggroup_VASPsol_issues_61-23issuecomment-2D1030676553&d=DwMFaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=CIoP38NMAXgBnnLsCELUobX_Pv7ewLJQLbNCmAz_GzOv1MUkwpdIa8DraR-2PBp4&s=72l4XrikgRx1amAg54H4aV1e4cYBPac0-EMqUTAfVQk&e=>, or unsubscribe<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_ACCOT2SB5S7TAEO7ZF6GSZ3UZVUTHANCNFSM5NSXVHOA&d=DwMFaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=CIoP38NMAXgBnnLsCELUobX_Pv7ewLJQLbNCmAz_GzOv1MUkwpdIa8DraR-2PBp4&s=Djmnu2vZ0Sp48arMB2pJ7159ToOCzzJ38e64BkKFOUI&e=>. Triage notifications on the go with GitHub Mobile for iOS<https://urldefense.proofpoint.com/v2/url?u=https-3A__apps.apple.com_app_apple-2Dstore_id1477376905-3Fct-3Dnotification-2Demail-26mt-3D8-26pt-3D524675&d=DwMFaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=CIoP38NMAXgBnnLsCELUobX_Pv7ewLJQLbNCmAz_GzOv1MUkwpdIa8DraR-2PBp4&s=gJIxaFkDht3xZgF43Uwvj4UJHU0Eok5aqAvXTkyhRdo&e=> or Android<https://urldefense.proofpoint.com/v2/url?u=https-3A__play.google.com_store_apps_details-3Fid-3Dcom.github.android-26referrer-3Dutm-5Fcampaign-253Dnotification-2Demail-2526utm-5Fmedium-253Demail-2526utm-5Fsource-253Dgithub&d=DwMFaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=CIoP38NMAXgBnnLsCELUobX_Pv7ewLJQLbNCmAz_GzOv1MUkwpdIa8DraR-2PBp4&s=tlnn43nRziba-XJbzdthNSrWMmEPwi_hW5SswJ4mkuk&e=>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[xnmao]

Dear Prof. @rhennig,

Sorry I missed them.
But I can't find a quadratic curve. The energy in OUTCAR is increased significantly by adding net charge.
Here is the full output:

	A	B	C	D	E
	net_charge	energy (in OUTCAR)	E-fermi (in OUTCAR)	FERMI_SHIFT (in stdout)	vacuum level (from LOCPOT)
0	-0.20	-93.781	-7.110	2.041	-2.031
1	-0.16	-93.494	-7.211	2.041	-2.033
2	-0.12	-93.204	-7.295	2.042	-2.035
3	-0.08	-92.910	-7.393	2.042	-2.037
4	-0.04	-92.612	-7.479	2.042	-2.039
5	0.00	-92.312	-7.562	2.041	-2.041
6	0.04	-92.007	-7.645	2.041	-2.044
7	0.08	-91.700	-7.718	2.042	-2.045
8	0.12	-91.389	-7.812	2.042	-2.047
9	0.16	-91.075	-7.894	2.041	-2.049
10	0.20	-90.757	-7.969	2.041	-2.051

Thank you for your kindness!

Best wishes,
Xinnan

[rhennig]

Dear Xinnan, Thanks for sharing the complete data. Looking at the results, it seems that the data you sent either already added the FERMI_SHIFT to the Fermi energy or the energy correction of FERMI_SHIFT times net_charge to the energy. The resulting quadratic curve looks correct. Best regards, Richard — Richard G. Hennig Professor and Alumni Professor of Materials Science and Engineering Associate Director, Quantum Theory Project Department of Materials Science and Engineering, University of Florida Email: ***@***.******@***.***> Web: http://hennig.mse.ufl.edu On Feb 5, 2022, at 4:35 PM, Xinnan Mao ***@***.******@***.***>> wrote: [External Email] Dear Prof. @rhennig<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_rhennig&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=eQ18oZzZ7NUJVSEGUTCEo3bhZgu5uvpV56UKSoLIezuw4H5PcGIVurUqXfjOxVLU&s=6y6f_6Vp2B43dcx4IDLSEaXQWl6IV72789jhd4uLeZ8&e=>, Sorry I missed them. But I can't find a quadratic curve. The energy in OUTCAR is increased significantly by adding net charge. Here is the full output: A B C D E net_charge energy (in OUTCAR) E-fermi (in OUTCAR) FERMI_SHIFT (in stdout) vacuum level (from LOCPOT) 0 -0.20 -93.781 -7.110 2.041 -2.031 1 -0.16 -93.494 -7.211 2.041 -2.033 2 -0.12 -93.204 -7.295 2.042 -2.035 3 -0.08 -92.910 -7.393 2.042 -2.037 4 -0.04 -92.612 -7.479 2.042 -2.039 5 0.00 -92.312 -7.562 2.041 -2.041 6 0.04 -92.007 -7.645 2.041 -2.044 7 0.08 -91.700 -7.718 2.042 -2.045 8 0.12 -91.389 -7.812 2.042 -2.047 9 0.16 -91.075 -7.894 2.041 -2.049 10 0.20 -90.757 -7.969 2.041 -2.051 Thank you for your kindness! Best wishes, Xinnan — Reply to this email directly, view it on GitHub<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_henniggroup_VASPsol_issues_61-23issuecomment-2D1030702777&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=eQ18oZzZ7NUJVSEGUTCEo3bhZgu5uvpV56UKSoLIezuw4H5PcGIVurUqXfjOxVLU&s=vxli0IuSugT0-bV25JODYWEUYzQ223fy9AlW0Cv1uT0&e=>, or unsubscribe<https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_ACCOT2Q4H52KNO56YLG4ZNDUZWJ2DANCNFSM5NSXVHOA&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=eQ18oZzZ7NUJVSEGUTCEo3bhZgu5uvpV56UKSoLIezuw4H5PcGIVurUqXfjOxVLU&s=tnKy6svlVkguO2WQ5jqNRYJ90gzI7gG7ZY4-675K44g&e=>. Triage notifications on the go with GitHub Mobile for iOS<https://urldefense.proofpoint.com/v2/url?u=https-3A__apps.apple.com_app_apple-2Dstore_id1477376905-3Fct-3Dnotification-2Demail-26mt-3D8-26pt-3D524675&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=eQ18oZzZ7NUJVSEGUTCEo3bhZgu5uvpV56UKSoLIezuw4H5PcGIVurUqXfjOxVLU&s=qojW2jxX8aNYusq0toZp5lvRkW75dP-L-Qk7XZOda4c&e=> or Android<https://urldefense.proofpoint.com/v2/url?u=https-3A__play.google.com_store_apps_details-3Fid-3Dcom.github.android-26referrer-3Dutm-5Fcampaign-253Dnotification-2Demail-2526utm-5Fmedium-253Demail-2526utm-5Fsource-253Dgithub&d=DwMCaQ&c=sJ6xIWYx-zLMB3EPkvcnVg&r=bSw49KO0hWxPjfCyM4jz2g&m=eQ18oZzZ7NUJVSEGUTCEo3bhZgu5uvpV56UKSoLIezuw4H5PcGIVurUqXfjOxVLU&s=Pdqjjz5hp0Dda_bma-I0q-f4saamKa3kp2VbBw0sHD0&e=>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[xnmao]

Dear @rhennig,

I may not have made it clear. I mean, if I plot column[B]-column[A]*column[D] against column[A], what I get seems to be a linear relationship between the grand-canonical electronic energy (?) and net charge as follows,

While Figure 2 in the VASPsol publication exhibits a typical quadratic behavior.

Best wishes,
Xinnan

[xnmao]

I have managed to do this now by applying a correction from the volume averaged electrostatic potential within the simulation cell, as described in equation 14, 15 of https://onlinelibrary.wiley.com/doi/10.1002/anie.200502540.

Many thanks!
Xinnan Mao

[jiangkaijiangkai]

Dear all,
I have a similar problem in repeating Figure 2 in the VASPsol publication (DOI:10.1063/1.5132354) A linear relationship between the energy and the electrode potential was obtained. After reading the paper (DOI:10.1002/anie.200502540) @xnmao recommended, I still could not get a quadratic relationship between the energy and the electrode potential after I added another correction. I guess it is about data processing. I've got to process the data as follows:

1. Electrode potential (U (V)):
   U=(E_vacuum_level(in LOCPOT) - E_Fermi(in OUTCAR) + FERMI_SHIFT(in stdout))/e.
   where e is 1. Similar to the xnamao'data above, the value of E_vacuum_level(in LOCPOT) is almost equal to the negative value of FERMI_SHIFT(in stdout)). So the U is actually equal to the negative value of E_Fermi (in OUTCAR).

2. The grand-canonical electronic energy (F (eV)):
   F=E_DFT(in OUTCAR) - n*U
   where n is the net charge, i.e., when the charge is added 0.02 to the neutral slab, the n = +0.02; when the charge is removed from the neutral system, the n = -0.02.

3. In the paper (DOI:10.1002/anie.200502540) @xnmao recommended, there is another correction is about the slab-background interaction: \int_{0}^{q} \left \langle Vtot\left ( Q \right ) \right \rangle dQ.
   where the Vtot(q) refers to the average potential in the unit cell.
   In my understanding, the Vtot(q) is the 'vacuum level' (in LOCPOT). However, I still could not get a quadratic relationship even after adding this correction. I also get confused about this correction, which I think it has been included in E_DFT(in OUTCAR) calculated in VASPsol.

Is there any mistake in my understanding? Could you give me more details?
Thanks!
Kai

[liam3935]

Dear Prof. @rhennig,

I'm wondering if geometry optimization is needed for each of the calculation set at different potential by NELECT? If so, will the quadratic behavior still hold if structure changes significantly?

Thanks,
Liam

[jiangkaijiangkai]

Hi @liam3935 ,
You can get the answer in #issue60 . And according to my limited experience, geometry optimization need to be done, and if structure changes significantly, you can NOT get the expected quadratic curve.
Kai

[rhennig]

The quadratic behavior should show if the structure is kept fixed. If the structure is allowed to relax, this will lead to an additional energy contribution from the structural changes that will obfuscate the quadratic dependence of the free energy on the charge. For infinitesimal charges, the change in structure should be small unless the system is near an instability. Have you calculated the free energy as a function of charge, making sure that you include the Q x V term? Richard From: Jiang Kai ***@***.***> Date: Thursday, April 20, 2023 at 8:55 PM To: henniggroup/VASPsol ***@***.***> Cc: Hennig, Richard ***@***.***>, Mention ***@***.***> Subject: Re: [henniggroup/VASPsol] calculation of the grand-canonical electronic energy (Issue #61) [External Email] Hi @liam3935<https://github.com/liam3935> , You can get the answer in #issue60 . And according to my limited experience, geometry optimization need to be done, and if structure changes significantly, you can NOT get the expected quadratic curve. Kai — Reply to this email directly, view it on GitHub<#61 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACCOT2RIEPTBJO5MJWMSPKLXCHLHRANCNFSM5NSXVHOA>. You are receiving this because you were mentioned.Message ID: ***@***.***>