UT theoretical chemistry code forum

I'm trying to analyze my CUBE files generated by SIESTA and have very strange behavior.

I've tested your H2O example it works fine.
My H2O file gave expected 0 on H and 8 on O.

Next I'm trying with simple thiol molecule HS-(CH2)4-CH3.
I've saved only partial core charges (used for NLCC) to make it as simple as possible i.e. very localized on atoms and exactly 0 between atoms and on the boundaries of the cell.

The resulting ACF file for 4A x 4A x 12A cell (molecule in the centre) with 98 x 98 x 172 grid:
# X Y Z CHARGE MIN DIST
----------------------------------------------------------------
1 0.0000 -0.1666 15.7422 0.0000 H 68.6933
2 0.0000 -1.1694 13.9206 0.0002 C 96.6529
3 1.6782 -2.4006 13.8791 0.0000 H 16.8121
4 -1.6782 -2.4006 13.8791 0.0000 H 17.3537
5 0.0000 0.6736 11.6876 0.0002 C 93.7153
6 1.6653 1.9249 11.8082 0.0000 H 22.1579
7 -1.6653 1.9249 11.8082 0.0000 H 22.6955
8 0.0000 -0.6715 9.1180 1.1506 C 96.7462
9 -1.6671 -1.9242 8.9995 0.0000 H 17.8655
10 1.6671 -1.9242 8.9995 0.0000 H 17.3273
11 0.0000 1.1447 6.8567 1.3867 C 94.2571
12 1.6697 2.3918 6.9525 0.0000 H 11.4817
13 -1.6697 2.3918 6.9525 0.0000 H 12.0206
14 0.0000 -0.2538 4.3246 1.3917 C 82.9846
15 -1.6876 -1.4600 4.1631 0.0000 H 6.5203
16 1.6876 -1.4601 4.1632 0.0000 H 5.8202
17 0.0001 2.0865 1.7261 5.4811 S 2.3419
18 0.0000 0.3746 -0.1867 0.0000 H 1.6147
----------------------------------------------------------------
NUMBER OF ELECTRONS: 9.41056

Don't be frightened with the non-integer charges - these are only partial cores.

Charge on H atoms is 0 as expected, charges on the 2 last C atoms and S atom are correct, while on the first C atoms it is practically 0, while they should be equal on all C atoms.

I've checked the CUBE file with Molekel and XCrysDen - the charges are exactly as I expect them.

BCF shows maxima as expected on S and C atoms:
# X Y Z CHARGE ATOM DISTANCE
-------------------------------------------------------------------------
1 0.0000 2.0054 1.7139 5.4811 17 0.0820
2 0.0000 -1.2341 13.8653 0.0002 2 0.0851
3 0.0000 -0.6171 9.0971 1.1506 8 0.0583
4 0.0000 -0.3085 4.3288 1.3917 14 0.0549
5 0.0000 0.6170 11.7119 0.0002 5 0.0616
6 0.0000 1.0798 6.9436 1.3867 11 0.1085
-------------------------------------------------------------------------
But shows 7 maxima instead of 6 (1S + 5C).

Even more strange are the min distances in ACF file. They are bigger than cell size and are growing in the direction of incorrectly analyzed atoms.

Grid refinement doesn't change the charges.
Input cell size change doesn't affect charges as well, but changes min distances in ACF (not in BCF).

ACF For 10A x 10A x 18A cell with same grid density:

# X Y Z CHARGE MIN DIST
----------------------------------------------------------------
1 0.0000 -0.1666 15.7422 0.0000 H 166.6554
2 0.0000 -1.1694 13.9206 0.0001 C 199.6380
3 1.6782 -2.4006 13.8791 0.0000 H 50.4000
4 -1.6782 -2.4006 13.8791 0.0000 H 49.7300
5 0.0000 0.6736 11.6876 0.0001 C 198.9569
6 1.6653 1.9249 11.8082 0.0000 H 69.4512
7 -1.6653 1.9249 11.8082 0.0000 H 68.7861
8 0.0000 -0.6715 9.1180 1.2040 C 201.7094
9 -1.6671 -1.9242 8.9995 0.0000 H 51.1873
10 1.6671 -1.9242 8.9995 0.0000 H 51.8531
11 0.0000 1.1447 6.8567 1.3939 C 199.2812
12 1.6697 2.3918 6.9525 0.0000 H 25.7107
13 -1.6697 2.3918 6.9525 0.0000 H 25.0439
14 0.0000 -0.2538 4.3246 1.4700 C 183.8613
15 -1.6876 -1.4600 4.1631 0.0000 H 7.5192
16 1.6876 -1.4601 4.1632 0.0000 H 8.5121
17 0.0001 2.0865 1.7261 5.4747 S 3.7647
18 0.0000 0.3746 -0.1867 0.0000 H 1.1343
----------------------------------------------------------------
NUMBER OF ELECTRONS: 9.54291


I've tried to visualize Bader volumes - but they are absolutely nonsensical - they don't have any symmetry of the system, instead each of the 7 volumes look like 9 ellipsoids in 3x3 regular pattern in one plane (as if it was square and then stretched to the shape of the cell), moreover all these 9 points are in one plane.

I haven't tried to draw Bader volumes for you H2O example since it is too big.

Analysis of just valence electrons and core+valence shows same behavior.

It seems like something is going wrong. Maybe we're not reading the file correctly. If you can send us the file, we'll see if something is going wrong with our cube reader and the seista files.

I've calculated your system with vasp, using the same grid density as yours. I translate your molecule a little, but the structure is the same as yours. I didn't have the same problem as you. My result is:

# X Y Z CHARGE MIN DIST Atom type
-------------------------------------------------------------------------------------------
1 7.5589 7.3923 19.7422 0.9771 0.8967 H
2 9.2371 5.1583 17.8791 0.9869 0.9350 H
3 5.8807 5.1583 17.8791 0.9873 0.9325 H
4 9.2242 9.4837 15.8081 0.9988 0.7109 H
5 5.8936 9.4837 15.8081 0.9991 0.9279 H
6 5.8918 5.6347 12.9995 0.9985 0.9278 H
7 9.2260 5.6347 12.9995 0.9989 0.9325 H
8 9.2286 9.9507 10.9525 0.9983 0.7130 H
9 5.8892 9.9507 10.9525 0.9985 0.9307 H
10 5.8713 6.0989 8.1631 0.9932 0.9438 H
11 9.2465 6.0988 8.1632 0.9937 0.9287 H
12 7.5589 7.9335 3.8133 0.9273 0.9413 H
13 7.5589 6.3895 17.9206 4.0775 1.0384 C
14 7.5589 8.2325 15.6876 4.0123 1.0874 C
15 7.5589 6.8874 13.1180 4.0183 1.0987 C
16 7.5589 8.7035 10.8567 4.0141 1.0853 C
17 7.5589 7.3051 8.3246 3.9880 1.0271 C
18 7.5589 9.6454 5.7261 6.0322 1.7583 S
NUMBER OF ELECTRONS: 38.00000

As you see, the result seems reasonable. Both H and C atoms have charges. I also tried the calculation without core charge. I could still get this result.

I think there's a problem with the charge density file you generate. For your cube file, it has 18 atoms but only has 6 significant maxima. That why some of the atoms don't have charges.

For the minimum distance problem, we've fixed that. Now the min dist seems reasonable.

Should I download a new version?

The wrong distances was not the main problem - I wouldn't care about them if the charges were correct.

Wenjie
- I've checked (visualized in 2 different programs) the densities before trying to get Bader charges - tehy are correct.
- there are six atoms with charges - (S and 5 C) because I was talking about core charges. Nevertheless, the program finds charges only on 3 or 4 atoms (depending on the cell size, which is very strange)
- if I take just valence electrons, I observe same strange behavior - (some C atoms stay without charge, and the total charge is incorrect as well). However, I see charge maxima on H atoms in visualization programs (since this is not O-H bond, but C-H)
- the results that you are showing - are they generated from your CHGCAR file? Have you tried to work with CUBE?
- Have you got my CUBE file? Did you try to visualize it? What results do you have with it? (I SHOULD have 6 significant maxima, but I get 7 reported by the program while actually only 3 or 4 atoms have nonzero charges)

I've checked the newer version (only main.f90 is newer - from March 8th).

I still get same wrong distances.
I also see that the bigger cell that I take, the more and more charge appears on C atoms (for 1.5x1.5A cell I see charge only on S, for 3x3A I see charge on S and three C atoms, for even bigger cells nothing changes anymore)
but the wronger distances as well.

By the way, I cannot work with more than 50x50x100 grid with double precision, so I need to switch to single precision. Although, I've checked of smaller grid sampling that it gives absolutely same results as dp.

We're still looking into the cube file problem. I think you have clearly shown that there is a problem interpreting the SEISTA->CUBE type file. We'll post back when we know what's going on, and update the code at the same time (with the min-dist fix included).

Thanks, Graeme.
Well, I wrote the export to CUBE myself, so I could make some errors as well.
But at least I can visualize it correctly in independent programs.
I'll try to look into it if you believe the error is there.

Ok, there was a silly mistake in the cube reader which was causing problems for non-cubic cells. I've updated the tar.gz file. The results look much better for the pentanethiol partial charges:

# X Y Z CHARGE MIN DIST
----------------------------------------------------------------
1 0.0000 -0.1666 15.7422 0.0000 8.2881
2 0.0000 -1.1694 13.9206 1.3877 9.8312
3 1.6782 -2.4006 13.8791 0.0000 4.1003
4 -1.6782 -2.4006 13.8791 0.0000 4.1658
5 0.0000 0.6736 11.6876 1.3708 9.6807
6 1.6653 1.9249 11.8082 0.0000 4.7072
7 -1.6653 1.9249 11.8082 0.0000 4.7640
8 0.0000 -0.6715 9.1180 1.3537 9.8360
9 -1.6671 -1.9242 8.9995 0.0000 4.2268
10 1.6671 -1.9242 8.9995 0.0000 4.1626
11 0.0000 1.1447 6.8567 1.3867 9.7086
12 1.6697 2.3918 6.9525 0.0000 3.3885
13 -1.6697 2.3918 6.9525 0.0000 3.4671
14 0.0000 -0.2538 4.3246 1.3917 9.1096
15 -1.6876 -1.4600 4.1631 0.0000 2.5535
16 1.6876 -1.4601 4.1632 0.0000 2.4125
17 0.0001 2.0865 1.7261 5.4812 1.5303
18 0.0000 0.3746 -0.1867 0.0000 1.2707
----------------------------------------------------------------
NUMBER OF ELECTRONS: 12.37181

Graeme, it work fine now, but it seems to me there is still some error - file with valence charge is not calculated properly, i.e. even total amount of electrons (the first and last H atoms get about 6 electrons), however, if I use total=core+valence it works correctly.
I'm sending you the files.

Also, after your last modifications I cannot read big files (even with single precision for q2), while in previous version this same files were opened correctly.
The problem is not even during calculations, but when it just tries to read it:

OPEN ... filename.cub
GAUSSIAN-STYLE INPUT FILE
Segmentation fault

Maybe you use memory somehow inefficiently? In my program of conversion of SIESTA to cube it stores in memory function on even much denser grids without any problems.

Yes, it was this segfault problem that prompted me to make the change which introduced the cube file bug. If you look in the cube_mod.f90, you try making the following change to the code. It just replaces an array multiplication with an explicit loop. For some reason, the array function has given me segfaults (but not all the time), whereas the explicit loop seems fine. Maybe the array function is creating a temporary variable?

-----------------
original code
-----------------

! GH: for some reason this is not working; replace with loop
chg%rho=chg%rho*vol
! DO n1=1,chg%npts(1)
! DO n2=1,chg%npts(2)
! DO n3=1,chg%npts(3)
! chg%rho(n1,n2,n3)=vol*chg%rho(n1,n2,n3)
! END DO
! END DO
! END DO

-------------------
replace with
-------------------

! GH: for some reason this is not working; replace with loop
! chg%rho=chg%rho*vol
DO n1=1,chg%npts(1)
DO n2=1,chg%npts(2)
DO n3=1,chg%npts(3)
chg%rho(n1,n2,n3)=vol*chg%rho(n1,n2,n3)
END DO
END DO
END DO