TURBOMOLE Modules > Molecular Properties, Wavefunction Analysis, and Interfaces to Visualization Tools

Smaller grid for calculation of electron density

(1/2) > >>

prasanta13:
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: ---$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

--- End code ---

The td.xyz starts as.

--- Code: ---#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
--- End code ---

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:
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: ---$pointval
grid1 points 200
grid2 points 200
grid3 points 200
--- End code ---

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:
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:
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?

uwe:
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.

Navigation

[0] Message Index

[#] Next page