integration results differ within twice run

[Ceasea]

Hi @aoterodelaroza

I run "yt nnm" twice in one Critic2 environment.
I found integration results of the two "yt nnm" are different.

Working flow is

Critic2 # enter command line mode
crystal CHGCAR
load CHGCAR
yt nnm # first run
yt nnm # second run

The two "yt nnm" provided different integration results.
(Maybe it's weird to run twice "yt nnm", but I did.)
'>.<'

My CHGCAR

CHGCAR.tar.gz

My compilier

GNU Fortran (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0

My compiling information

critic2 (development), commit (config) no-git
compile host: Linux ceasea-VirtualBox 4.15.0-47-generic #50-Ubuntu SMP Wed Mar 13 10:40:37 UTC 2019 i686 i686 i686 GNU/Linux
compile date: 2019年 07月 11日 星期四 16:32:59 CST
using f77: gfortran -g -O2
f90: gfortran -g -O2 -ffree-form -ffree-line-length-none -fopenmp
c: gcc -g -O2 -Iqhull/
ldflags:
debug?: no
compiled dat: /usr/local/share/critic2
datadir: /usr/local/share/critic2
dic file: /usr/local/share/critic2/cif/cif_core.dic
...was found?: T

[aoterodelaroza]

Thanks for the report. I think I know why that happens. In your crystal, you have NNMs (which could originate spuriously from using the pseudo-density to calculate the basins). When you run YT for the first time, those NNMs are added to the critical point list for the field as maxima, and therefore act as "atoms" in the second run.

In general, whenever YT finds a new attractor it is assigned to the nearest atom using a distance criterion (the RATOM keyword (1 bohr by default). This is done to prevent noisy densities from spawning great very many maxima that would both slow down the calculation and make the output a lot more confusing. In addition, every maximum in the critical point list of the field is listed in the YT output, even if it integrated to zero. This is done to avoid "losing" atoms in the process.

In your case, the difference between the two runs arises from the fact that, in the first run, your list of maxima contains only the Al atom whereas in the second run it contains the Al and all the maxima found in the first run. Since you have so many maxima and the system is so small, I would suggest using:

YT NNM NOATOMS

In this way, you will get the raw output from YT and there will be no accumulation of basin properties based on distance criteria. If you do that, you'll see that the run gives the same result both times.

I'm not sure if this is a bug, since YT is doing what it was designed to do. However, I do agree that it is unsettling that two runs on the same thing give different results.

[Ceasea]

! Pre-allocate atoms as maxima
allocate(bas%xattr(3,s%c%ncel))
bas%xattr = 0d0
if (bas%atexist) then
bas%nattr = s%c%ncel
do i = 1, s%c%ncel
bas%xattr(:,i) = s%c%atcel(i)%x
end do
else
bas%nattr = 0
end if

This code is from yt@proc.f90 from line 63 to line 73.

bas%xattr, which stores the coordinations of maxima, and bas%nattr, which stores the number of
maxima, are both initialized at the beginning of yt proc. ( if I was wrong, correct me)

I was wondering if the result of the first run would influence in the second one since bas%xattr and bas%nattr are initialized, or I was missing some important code in other files?

[aoterodelaroza]

That is correct. Both bas%nattr and bas%xattr are initialized to the number of atoms. And the maxima of YT are saved to the field CP list in integration@proc.f90, routine int_reorder_gridout. But right now I don't see right now the code where it recieves the CPs from the reference field. Hmmm. I'll take a look at the end of today and come back to you.

[aoterodelaroza]

This is the part of the code that does the deed (in integration@proc.f90, int_reorder_gridout):

if (bas%atexist) then
    do i = 1, nattr0
       nid = ff%c%identify_atom(bas%xattr(:,i),icrd_crys,distmax=bas%ratom)
       if (nid > 0) then
          assigned(i) = nid
       else
          ! maybe the closest point is a known nnm
          call ff%nearest_cp(bas%xattr(:,i),nid,dist,type=ff%typnuc)
          if (dist < bas%ratom) then
             assigned(i) = nid
          end if
       end if
    end do
end if

In the second pass, maxima can still not be assigned to adjacent atoms (because ratom = 1bohr and you only have Al) but they can be assigned to the NNMs found in the first run, and so they are.

[aoterodelaroza]

The bas%atexists can be made .false. by using the NOATOM keyword, as mentioned above.

[Ceasea]

Hi @aoterodelaroza,

Combining your explanation, I read the code and found the problem, and solved the problem basically in my tests.
I have pulled a request.

[aoterodelaroza]

Hi @Ceasea, Sorry for the wait - I have been all week at a conference, and will be until Monday. I do not think this is such idea because it makes the field inconsistent. The CP list should always be populated by at least the atoms, in the same order as the crystal structure. In addition, there are two lists that are coordinated: one with and one without symmetry. Finally, this fix also breaks the encapsulation - basically that if you changing the CP list anywhere there will come a time when it will be impossible to keep track of who touches the CP list and why. I have just pushed a change to the development version where, instead of adding the list of atoms to the initial YT and BADER lists of attractors, it adds all the NNMs in the field. This should make the two runs consistent with each other, even though the nnm names in the first list will be "??" and something like "n230" in the second. Can you please git pull and check?

[Ceasea]

yes，after i pulled the request, i realized that my changes are too violent and may cause other problems that i donot know. i found that in command-line mode,when i excuted unload command  between two YT NNM run, the results would be the same. However, repeating loading fields will cost much time if these field are big and is not beautiful. So i checked the unload module and found that parameter works only for this problem. Anyway, i will try the newest version. thanks. 汪进凯 邮箱：ceasea@i.shu.edu.cn 签名由 网易邮箱大师 定制 On 07/27/2019 22:20, aoterodelaroza wrote: Hi @Ceasea, Sorry for the wait - I have been all week at a conference, and will be until Monday. I do not think this is such idea because it makes the field inconsistent. The CP list should always be populated by at least the atoms, in the same order as the crystal structure. In addition, there are two lists that are coordinated: one with and one without symmetry. Finally, this fix also breaks the encapsulation - basically that if you changing the CP list anywhere there will come a time when it will be impossible to keep track of who touches the CP list and why. I have just pushed a change to the development version where, instead of adding the list of atoms to the initial YT and BADER lists of attractors, it adds all the NNMs in the field. This should make the two runs consistent with each other, even though the nnm names in the first list will be "??" and something like "n230" in the second. Can you please git pull and check? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[Ceasea]

Hi @aoterodelaroza
I have checked new version.
The problem is solved. Thanks.