Ground state wavefunction of H atom not spherical symmetric?

Total energy, geometry optimization, DFT+U, spin....

Moderator: bguster

Locked
kingkingred
Posts: 2
Joined: Wed Feb 05, 2014 1:47 am

Ground state wavefunction of H atom not spherical symmetric?

Post by kingkingred » Wed Feb 05, 2014 2:10 am

Dear all,
I just copied the code in tutorial 2.4 (http://www.abinit.org/documentation/helpfiles/for-v7.6/tutorial/lesson_base2.html#24), in which dataset1 deals with H2 molecule and dataset2 deals with an isolated H atom. I used xcrysden to plot the real part of WFK file of dataset2. It's completely not spherically symmetric. It's hard to understand since the nondegenerate 1s ground state (even for the real part) should be spherically symmetric. However, other results seem to be good. And the density is indeed spherical symmetric. Could you please help me on this?

Code: Select all

# H2 molecule in a big box
#
# This file to optimize the H2 bond length, compute the associated total
# energy, then to compute the total energy of the isolated H atom.
# Here, the ecut and acell are fixed : the double loop reduces
# effectively to a single loop.

ndtset 2  udtset 1 2

#Definition of the unit cell and ecut,
#for which one will have to make a convergence study
 acell:? 12 12 12  acell+? 2 2 2
 ecut 30

#First dataset : find the optimal bond length of H2, and associated total energy
   natom?1 2             # There are two atoms
  ionmov?1 3             # Use the modified Broyden algorithm
   ntime?1 10            # Maximum number of Broyden "timesteps"
  tolmxf?1 5.0d-4        # Stopping criterion for the geometry optimization : when
                         # the residual forces are less than tolmxf, the Broyden
                         # algorithm can stop
   xcart?1 -0.7  0.0 0.0 # The starting values of the
            0.7  0.0 0.0 # atomic coordinates
  toldff?1 5.0d-5        # Will stop the SCF cycle when, twice in a row,
                         # the difference between two consecutive evaluations of
                         # forces differ by less than toldff (in Hartree/Bohr)
   nband?1  1            # Just one band

#Second dataset : get the total energy of the isolated atom
   natom?2 1             # There is one atom
  nsppol?2 2             # Spin-polarized calculation
  occopt?2 2             # Allow occupation numbers to be set by hand
   nband?2 1 1           # Number of bands for spin up and spin down
     occ?2 1.0  0.0      # Occupation numbers for spin up state and spin down state.
  toldfe?2 1.0d-6        # Will stop the SCF cycles when, twice in a row,
                         # the difference between two consecutive evaluations
                         # of total energy differ by less than toldfe (in Hartree)
                         # This value is way too large for most realistic studies of materials
   xcart?2 0.0 0.0 0.0   # The atom is located at the origin
  spinat?2 0.0 0.0 1.0   # Initialisation of spin

#rprim 1 0 0  0 1 0  0 0 1 # This line, defining orthogonal primitive vectors,
                           # is commented, because it is precisely the default value of rprim

#Definition of the atom types
ntypat 1          # There is only one type of atom
znucl 1           # The keyword "znucl" refers to the atomic number of the
                  # possible type(s) of atom. The pseudopotential(s)
                  # mentioned in the "files" file must correspond
                  # to the type(s) of atom. Here, the only type is Hydrogen.
                         

#Definition of the atoms
typat 1 1         # For the first dataset, both numbers will be read,
                  # while for the second dataset, only one number will be read

#Definition of the k-point grid
kptopt 0          # Enter the k points manually
nkpt 1            # Only one k point is needed for isolated system,
                  # taken by default to be 0.0 0.0 0.0

#Definition of the SCF procedure
nstep 10          # Maximal number of SCF cycles
#toldfe is no more defined, as toldff is used above...
diemac 1.0        # Although this is not mandatory, it is worth to
                  # precondition the SCF cycle. The model dielectric
                  # function used as the standard preconditioner
                  # is described in the "dielng" input variable section.
                  # Here, we follow the prescriptions for molecules
                  # in a big box
# add to conserve old < 6.7.2 behavior for calculating forces at each SCF step
 optforces 1

User avatar
gmatteo
Posts: 291
Joined: Sun Aug 16, 2009 5:40 pm

Re: Ground state wavefunction of H atom not spherical symmet

Post by gmatteo » Thu Feb 06, 2014 2:33 pm

I used xcrysden to plot the real part of WFK file of dataset2.


How did you extract the real part?
I used cut3d and vesta to plot the norm of the 1s state and it's spherically symmetric as expected (see attached file)
Besides, in dataset12, one has n(r) = |psi_{1s}(r)|^2 so it seems strange that you have a spherically symmetric density whereas psi_{1s} is not!
Attachments
Screen shot 2014-02-06 at 14.29.50.png

kingkingred
Posts: 2
Joined: Wed Feb 05, 2014 1:47 am

Re: Ground state wavefunction of H atom not spherical symmet

Post by kingkingred » Thu Feb 06, 2014 2:50 pm

gmatteo wrote:
I used xcrysden to plot the real part of WFK file of dataset2.


How did you extract the real part?
I used cut3d and vesta to plot the norm of the 1s state and it's spherically symmetric as expected (see attached file)
Besides, in dataset12, one has n(r) = |psi_{1s}(r)|^2 so it seems strange that you have a spherically symmetric density whereas psi_{1s} is not!



cut3d just provides the 'real part' option when it processes the unformatted WFK file.
Indeed, that's why I posted it here.

User avatar
gmatteo
Posts: 291
Joined: Sun Aug 16, 2009 5:40 pm

Re: Ground state wavefunction of H atom not spherical symmet

Post by gmatteo » Thu Feb 06, 2014 3:33 pm

There was a bug if that particular option was selected (Xcrysden uses general grids, see http://www.xcrysden.org/doc/XSF.html#scalar-field)

This patch should solve the problem. Alternatively, you can get the new version of wffile.F90 attached to this reply.
(I had to replace the F90 extension with log to upload the file. Replace .log with .F90, put the file in src/83_cut3d and recompile with make)

Code: Select all


=== modified file 'src/83_cut3d/wffile.F90'
--- src/83_cut3d/wffile.F90     2013-01-23 13:57:47 +0000
+++ src/83_cut3d/wffile.F90     2014-02-06 14:15:20 +0000
@@ -950,6 +950,7 @@
          shift_tau(:) = 0.0
          read (*,*) outputchar
          if (outputchar == 'y' .or. outputchar == 'Y') then
+           MSG_ERROR("Shift is buggy, don't use it")
            write(std_out,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,') :'
            write(std_out,*)
            read (*,*) gridshift1, gridshift2, gridshift3
@@ -1109,6 +1110,7 @@
          shift_tau(:) = 0.0
          read (*,*) outputchar
          if (outputchar == 'y' .or. outputchar == 'Y') then
+           MSG_ERROR("Shift is buggy, don't use it")
            write(std_out,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,') :'
            write(std_out,*)
            read (*,*) gridshift1, gridshift2, gridshift3
@@ -1165,16 +1167,19 @@
          write(unout,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D'
          write(unout,*) 'datagrids'
          write(unout,'(1X,A)') 'DATAGRID_3D_DENSITY'
-         write(unout,*) nr1,nr2,nr3
+         write(unout,*) nr1+1,nr2+1,nr3+1
          write(unout,*) '0.0 0.0 0.0 '
          do ir1 = 1,3
            write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3)
          end do
 
-         do ir3=1,nr3
-           do ir2=1,nr2
-             do ir1=1,nr1
-               write(unout,'(ES17.10)') fofr(1,ir1,ir2,ir3)
+         do ir3=1,nr3+1
+           ii3=mod(ir3-1,nr3) + 1
+           do ir2=1,nr2+1
+             ii2=mod(ir2-1,nr2) + 1
+             do ir1=1,nr1+1
+               ii1=mod(ir1-1,nr1) + 1
+               write(unout,'(ES17.10)') fofr(1,ii1,ii2,ii3)
              end do
            end do
          end do
Attachments
wffile.log
(46.1 KiB) Downloaded 240 times

Locked