Author Topic: Smaller grid for calculation of electron density  (Read 192 times)

prasanta13

  • Full Member
  • ***
  • Posts: 34
  • Karma: +0/-0
Smaller grid for calculation of electron density
« on: July 04, 2022, 10:34:53 AM »
Hi there, I am trying to calculate the electron density of H2 molecule. I have ran a simple scf calculation. Control file is provided below. I got the output as td.xyz.
Code: [Select]
$title
$symmetry c1
$user-defined bonds    file=coord
$coord    file=coord
$optimize
 internal   off
 redundant  off
 cartesian  on
 global     off
 basis      off
$atoms
h  1-2                                                                         \
   basis =h def2-QZVPPD
$basis    file=basis
$scfmo   file=mos
$closed shells
 a       1                                      ( 2 )
$scfiterlimit       30
$scfconv        7
$thize     0.10000000E-04
$thime        5
$scfdamp   start=0.300  step=0.050  min=0.100
$scfdump
$scfintunit
 unit=30       size=0        file=twoint
$scfdiis
$maxcor    500 MiB  per_core
$scforbitalshift  automatic=.1
$drvopt
   cartesian  on
   basis      off
   global     off
   hessian    on
   dipole     on
   nuclear polarizability
$interconversion  off
   qconv=1.d-7
   maxiter=25
$coordinateupdate
   dqmax=0.3
   interpolate  on
   statistics    5
$forceupdate
   ahlrichs numgeo=0  mingeo=3 maxgeo=4 modus=<g|dq> dynamic fail=0.3
   threig=0.005  reseig=0.005  thrbig=3.0  scale=1.00  damping=0.0
$forceinit on
   diag=default
$energy    file=energy
$grad    file=gradient
$forceapprox    file=forceapprox
$pointval fmt=xyz
$rundimensions
   natoms=2
   nbf(CAO)=76
   nbf(AO)=66
$last step     dscf
$orbital_max_rnorm 0.44543416007246E-05
$last SCF energy change = -1.1335538
$charge from dscf
          0.000 (not to be modified here)
$dipole from dscf
  x     0.00000000000002    y    -0.00000000000017    z     0.00000000000088    a.u.
   | dipole | =    0.0000000000  debye
$end

The td.xyz starts as.
Code: [Select]
#origin           0.000000      0.000000      0.000000
#vector1          1.000000      0.000000      0.000000
#vector2          0.000000      1.000000      0.000000
#vector3          0.000000      0.000000      1.000000
#grid1  start   -5.313977  delta    0.205059  points     44
#grid2  start   -7.744278  delta    0.203075  points     46
#grid3  start    4.056177  delta    0.205458  points     40
#title for this grid 101
#density
#plotdata
# cartesian coordinates x,y,z and f(x,y,z)
       -5.31397670         -7.74427833          4.05617653          0.00000006
       -5.10891784         -7.74427833          4.05617653          0.00000008
       -4.90385899         -7.74427833          4.05617653          0.00000010
       -4.69880013         -7.74427833          4.05617653          0.00000013
       -4.49374127         -7.74427833          4.05617653          0.00000017
       -4.28868241         -7.74427833          4.05617653          0.00000021
       -4.08362355         -7.74427833          4.05617653          0.00000026
       -3.87856469         -7.74427833          4.05617653          0.00000032
       -3.67350583         -7.74427833          4.05617653          0.00000039

What I get is that delta is the grid spacing, for each direction, x, y and z. Can I obtain a smaller grid spacing? Rather than 0.20  I want to get 0.02 or even smaller?

antti_karttunen

  • Sr. Member
  • ****
  • Posts: 221
  • Karma: +1/-0
Re: Smaller grid for calculation of electron density
« Reply #1 on: July 06, 2022, 09:41:36 AM »
Hi,

yes, you can set the density of the grid yourself. This is done with sub-keywords grid1, grid2, and grid3 under the $pointval keyword. The option "points" determines the number of grid points in each direction:

Code: [Select]
$pointval
grid1 points 200
grid2 points 200
grid3 points 200

You can also adjust the grid vectors and ranges this way. More information in the manual section 22.2.26 Keywords for wave function analysis and generation of plotting data.

Best,
Antti

uwe

  • Global Moderator
  • Hero Member
  • *****
  • Posts: 515
  • Karma: +0/-0
Re: Smaller grid for calculation of electron density
« Reply #2 on: July 07, 2022, 04:09:48 PM »
Hi,
for a smaller grid size and if you do not need a specific number, just add 'quick' to the $pointval line.
The number of points is set to half of the default in each direction. Might be too sparse, but TmoleX uses it when you generate 3D grid data for visualization by clicking on the 'quick' button instead of the 'high-res' one. Usually sufficient for visualization purposes.
Best, Uwe

prasanta13

  • Full Member
  • ***
  • Posts: 34
  • Karma: +0/-0
Re: Smaller grid for calculation of electron density
« Reply #3 on: July 08, 2022, 07:51:14 AM »
Actually, not for plotting purposes but for numerical reasons, I am trying to integrate the total electron density so that it matches the total electron number of the molecule.
For H2 molecule, the total density value should be 2.
Without changing the grid start, and origin I am just increasing the grid points (thus decreasing the delta or grid size) such that I reach the value closest to 2. I am showing the results,
delta density
0.03 1.976
0.025 1.980
0.015 1.9869

Are there any quick hacks to get the electron density very close to the actual one? Nay special keyword to write in control?
« Last Edit: July 08, 2022, 02:33:36 PM by prasanta13 »

uwe

  • Global Moderator
  • Hero Member
  • *****
  • Posts: 515
  • Karma: +0/-0
Re: Smaller grid for calculation of electron density
« Reply #4 on: July 08, 2022, 11:23:01 AM »
Hi,

two things you could/might have to do:

1. grid

To catch all of the density you probably need to use a bigger box which reaches further out in space - depending on the basis set.

(a bit off topic):

To find out which atomic orbitals contribute to the density (by checking the occupied molecular orbitals and printing the contribution of each atomic orbital) try:

tm2molden mostat <list-of-MOs> above <threshold>

So if you have 40 occupied orbitals (eiger quickly tells you how many you have) and you want to see all AOs which contribute to each MO more than, say,  0.01% just call on the command line:

tm2molden mostat 1-40 above 0.01

Enter a few times <enter> and then you will get a file named mostat.out which looks like that:

 MO Index (from eiger):       3
 Sym=     3a
 Energy [Hartree]= -.99177169241775E+01
 Energy [eV]= -.26987482366088E+03
 Spin= Alpha
 Occup= 2.000000
contribution of AO   2 c     1 s        is:   24.939 %
contribution of AO   6 c     1 s        is:   24.939 %
contribution of AO   3 c     1 s        is:   24.939 %
contribution of AO   5 c     1 s        is:   24.939 %
contribution of AO   2 c     2 s        is:    0.052 %

and so on. That is useful to see how many diffuse AOs contribute to the occupied orbitals and thus to the density.

2. numerical accuracy

Another source of numerical inaccuracy is the format of either plt or cube files. They are by default single-precision data. As the density is written to a grid and then summed up, just from pure numerics the sum will never exactly be the number of electrons (well, never is not correct, it could coincidentally be just an integer number...).

Try to add fmt=dtx to the $pointval input line instead. This will print double precision data as ASCII text (dtx = double precision text file). Looks like that:

#origin           0.000000      0.000000      0.000000
#vector1          1.000000      0.000000      0.000000
#vector2          0.000000      1.000000      0.000000
#vector3          0.000000      0.000000      1.000000
#grid1  start   -5.767470  delta    0.202135  points     66
#grid2  start   -5.767470  delta    0.202829  points     59
#grid3  start   -5.020450  delta    0.204532  points     67
#title for this grid 101
#density
#plotdata
# cartesian coordinates x,y,z and f(x,y,z)
-5.7674700000000E+00 -5.7674700000000E+00 -5.0204500000000E+00  3.9793580979130E-10
-5.5653350769231E+00 -5.7674700000000E+00 -5.0204500000000E+00  5.8733532389287E-10
-5.3632001538462E+00 -5.7674700000000E+00 -5.0204500000000E+00  8.5021125663768E-10
-5.1610652307692E+00 -5.7674700000000E+00 -5.0204500000000E+00  1.2071193363237E-09
...

The last column then needs to be summed up. Whether or not that will increase accuracy - it's at least worth a try.

prasanta13

  • Full Member
  • ***
  • Posts: 34
  • Karma: +0/-0
Re: Smaller grid for calculation of electron density
« Reply #5 on: July 10, 2022, 04:01:53 PM »
Thanks for the info UWE,
First I tried the 2nd method prescribed by you but I do not see any significant change up to 3rd decimal.
Using 0.015 grid with single and double precision leads to the same density 1.987 for the H2 molecule.
I tried with a smaller grid 0.005, the output file td.xyz has become 382 GB, which I can not process.

Thanks again
Prasanta

prasanta13

  • Full Member
  • ***
  • Posts: 34
  • Karma: +0/-0
Re: Smaller grid for calculation of electron density
« Reply #6 on: July 11, 2022, 01:11:21 PM »
Hi there,
Again this issue I am finding difficult to digest.
I am now issuing the grid start points, delta and points, but the program does not utilize it.
Sharing the control file.
Code: [Select]
$title
$symmetry c1
$user-defined bonds    file=coord
$coord    file=coord
$optimize
 internal   off
 redundant  off
 cartesian  on
 global     off
 basis      off
$atoms
br 1-2                                                                         \
   basis =br def2-QZVPPD
$basis    file=basis
$scfmo   file=mos
$closed shells
 a       1-35                                   ( 2 )
$scfiterlimit       30
$scfconv        7
$thize     0.10000000E-04
$thime        5
$scfdamp   start=0.300  step=0.050  min=0.100
$scfdump
$scfintunit
 unit=30       size=0        file=twoint
$scfdiis
$maxcor    500 MiB  per_core
$scforbitalshift  automatic=.1
$drvopt
   cartesian  on
   basis      off
   global     off
   hessian    on
   dipole     on
   nuclear polarizability
$interconversion  off
   qconv=1.d-7
   maxiter=25
$coordinateupdate
   dqmax=0.3
   interpolate  on
   statistics    5
$forceupdate
   ahlrichs numgeo=0  mingeo=3 maxgeo=4 modus=<g|dq> dynamic fail=0.3
   threig=0.005  reseig=0.005  thrbig=3.0  scale=1.00  damping=0.0
$forceinit on
   diag=default
$energy    file=energy
$grad    file=gradient
$forceapprox    file=forceapprox
$pointval fmt=xyz
grid1 start -7.309453 delta 0.025 points 600
grid2 start -7.337762 delta 0.025 points 578
grid3 start -9.606856 delta 0.025 points 698
$rundimensions
   natoms=2
   nbf(CAO)=242
   nbf(AO)=196
$last step     dscf
$orbital_max_rnorm 0.14281617249944E-04
$last SCF energy change = -5144.9208
$charge from dscf
          0.000 (not to be modified here)
$dipole from dscf
  x    -0.00000000000081    y    -0.00000000000163    z     0.00000000002117    a.u.
   | dipole | =    0.0000000001  debye
$end
The header of the td.xyz shows this,
Code: [Select]
#origin           0.000000      0.000000      0.000000
#vector1          1.000000      0.000000      0.000000
#vector2          0.000000      1.000000      0.000000
#vector3          0.000000      0.000000      1.000000
#grid1  start   -5.241196  delta    0.017071  points    600
#grid2  start   -6.029632  delta    0.020265  points    578
#grid3  start   -1.012296  delta    0.012420  points    698
#title for this grid 101
#density
#plotdata
The points are kept fixed while delta and grid start points are different. Any help?
Thanks in advance
Prasanta