TURBOMOLE Users Forum
TURBOMOLE Modules => Molecular Properties, Wavefunction Analysis, and Interfaces to Visualization Tools => Topic started by: prasanta13 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.
$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.
#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?
-
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:
$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
-
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
-
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?
-
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.
-
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
-
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.$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,
#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