How to add a proton in the system

[nidhimangal]

Hi Francois,

I am a new Qbox user. Thank a lot for helping me. My Query is that if we include net_charge value that will change the net_charge value of the system. In my case I want to calculate the single proton surrounded by water molecules. How can it specifically change the charge of that particular hydrogen (proton). If I change the valance charge value (number of valence electron) from 1 to 0 in xml file of hydrogen species. Is that wrong procedure. Please guide me. Thank you

Kind regards
Nidhi

[fgygi]

Hi Nidhi,
(I insert here the first part of the answer and discuss it further below)

To add a proton to a system, you should add a hydrogen atom and then remove one electron from the system by changing the net_charge value:

Code: Select all

# add a proton
# define the hydrogen species
species hydrogen http://fpmd.ucdavis.edu/potentials/H/H_HSCV_PBE-1.0.xml
# add a hydrogen atom
atom H hydrogen  1.0 2.0 3.0
# remove an electron
set net_charge 1

For example, here is a script that will compute the ground state of an H2+ ion:

Code: Select all

# H2+ ion
set cell 20 0 0  0 20 0  0 0 20
species hydrogen http://fpmd.ucdavis.edu/potentials/H/H_HSCV_PBE-1.0.xml
atom H1 hydrogen 0.00 0.00 0.00
atom H2 hydrogen 2.36 0.00 0.00
set net_charge 1
set ecut 35
set xc PBE
set wf_dyn PSD
randomize_wf
run 0 50

Note that it is not possible to specify on which atom the charge will be located in a DFT calculation. This will be dictated by the solution of the Kohn-Sham equations. It is only possible to specify the total number of electrons in the calculation. In Qbox, the total number of electrons is calculated using the valence charge in the pseudopotential files, and it can be adjusted using the net_charge variable.

For the system you consider (water molecules + proton) it is very likely that the ground state electron density will be mostly located on the water molecules and that the additional hydrogen atom will be a bare proton. For example, in the case of the (H2O--H--OH2)+ cluster (total charge=16), it appears that the density is mostly on the two water molecules and the hydrogen in the middle has practically no charge (and it is therefore just a proton).

It is not possible to change the number of electrons by modifying the valence charge in the pseudopotential files. Changing the value of <valence_charge> in a pseudopotential file would yield incorrect results, since Qbox expects the potential given in the file to behave asymptotically as -<valence_charge>/r at large distance.

Here is an input file that works for the (H2O--H--OH2)+ ion. It should be easy to modify it to move the molecules, and/or add other molecules:

Code: Select all

# ground state of the (H2O--H--OH2)+ ion
set cell  20 0 0  0 20 0  0 0 20
species oxygen http://fpmd.ucdavis.edu/potentials/O/O_HSCV_PBE-1.0.xml
species hydrogen http://fpmd.ucdavis.edu/potentials/H/H_HSCV_PBE-1.0.xml
atom O1 oxygen   -2.22734253  -0.01179045   0.40254298
atom O2 oxygen    2.30366467   0.05714217  -0.39447360 
atom H1 hydrogen   -3.13951043   1.47354644  -0.22134861
atom H2 hydrogen   -3.10089639  -1.52074908  -0.21985775
atom H3 hydrogen    3.17936306   1.54998719   0.26291380
atom H4 hydrogen    3.21673943  -1.44331112   0.19061675
atom H5 hydrogen    0.03826565   0.01958876   0.00451500
set net_charge +1
set ecut 85
set wf_dyn PSDA
randomize_wf
set xc PBE
run 0 150
save gs1.xml

Note that in the above example, the geometry given in the input file is already a well optimized geometry. If you modify the position of the water molecules or add other molecules, you will have to reoptimize the geometry using e.g. CG steps.

Francois