basin.f90

!!------ Forming basin for specific real space function, and integrate real space functions in the basins
!!------ The method is adapted from J. Phys.: Condens. Matter 21 (2009) 084204
!Grid must be ortho
subroutine basinana
use defvar
use basinintmod
use util
use GUI
implicit real*8(a-h,o-z)
integer walltime1,walltime2
integer :: igridmethod=3
integer,allocatable :: attconvold(:),realattconv(:),usercluslist(:),tmparr(:)
character basinfilename*200,selectyn,c80tmp*80,ctmp1*20,ctmp2*20,ctmp3*20,ctmp4*20,c2000tmp*2000,c200tmp*200
real*8 :: threslowvalatt=1D-5
real*8,allocatable :: tmparrf(:)
ishowattlab=1
ishowatmlab=0
ishowatt=1
!Don't show topology analysis results to avoid confusion
ishowCPlab=0
ishowpathlab=1
ishow3n3=0
ishow3n1=0
ishow3p1=0
ishow3p3=0
!When enter this module, we reset the whole state, which is signified by numatt. Because the user may calculate grid data at external functions,
!and hence messed up grid setting (orgx/y/z,dx/y/z...), so all of the functions in this module that rely on grid setting will be totally wrong
numatt=0

call delvirorb(1)

do while(.true.)
	write(*,*)
	write(*,*) "                 ============= Basin analysis ============="
	write(*,*) "-10 Return to main menu"
	if (numatt>0) write(*,*) "-6 Set parameter for attractor clustering or manually perform clustering"
	if (numatt==0) write(*,*) "-6 Set parameter for attractor clustering"
	if (numatt>0) write(*,*) "-5 Export basins as cube file"
	if (numatt>0) write(*,*) "-4 Export attractors as pdb/pqr/txt file"
	if (numatt>0) write(*,*) "-3 Show information of attractors"
	if (numatt>0) write(*,*) "-2 Measure distances, angles and dihedral angles between attractors or atoms"
	if (igridmethod==1) write(*,*) "-1 Select the method for generating basins, current: On-grid"
	if (igridmethod==2) write(*,*) "-1 Select the method for generating basins, current: Near-grid"
	if (igridmethod==3) write(*,*) "-1 Select the method for generating basins, current: Near-grid with refinement"
	if (numatt>0) write(*,*) " 0 Visualize attractors and basins"
	if (numatt>0) then
		write(*,*) " 1 Regenerate basins and relocate attractors"
	else
		write(*,*) " 1 Generate basins and locate attractors"
	end if
	if (numatt>0) then
		write(*,*) " 2 Integrate real space functions in the basins"
		write(*,*) " 3 Calculate electric multipole moments in the basins"
		write(*,*) " 4 Calculate localization index and delocalization index for the basins"
		write(*,*) " 5 Output orbital overlap matrix in basins to BOM.txt in current folder"
! 		write(*,*) " 100 Integrate real space functions in the basins with multi-level refinement"
		if (ifuncbasin==1) then
		    write(*,*) " 6 Output orbital overlap matrix in atoms to AOM.txt in current folder"
			write(*,*) " 7 Integrate real space functions in AIM basins with mixed type of grids"
			write(*,*) " 8 Calculate electric multipole moments in AIM basins with mixed type of grids"
			write(*,*) " 9 Obtain atomic contribution to population of external basins"
        else if (ifuncbasin==9) then
            write(*,"(a)") " 10 Calculate high ELF localization domain population and volume (HELP, HELV)"
		end if
		write(*,*) "11 Calculate orbital compositions contributed by various basins"
	end if
	read(*,*) isel
	if (numatt==0.and.(isel==-5.or.isel==-4.or.isel==-3.or.isel==-2.or.isel==0.or.isel==2.or.isel==3.or.isel==4.or.isel==5)) then
		write(*,*) "Error: You should use option 1 to generate basins first!"
		cycle
	end if

	if (isel==-10) then
		idrawbasinidx=-10 !Don't display interbasin in any other GUI
		return
		
	else if (isel==-6) then
		do while(.true.)
			write(*,*) "0 Return"
			write(*,"(a,f14.5,a)") " 1 Set relative value difference criterion, current:",valcritclus*100,"%"
			write(*,"(a,i3)") " 2 Set the multiplier for distance criterion, current:",mergeattdist
			if (numatt>=2) write(*,*) "3 Cluster specified attractors"
			read(*,*) isel2
			if (isel2==0) then
				exit
			else if (isel2==1) then
				write(*,"(a)") " Input a value in percentage, e.g. 0.5"
				write(*,"(a)") " Note: Inputting X means if the value difference between two attractors is less than X%, &
				then they will be regarded as degenerate and may be clustered according to distance criterion"
				read(*,*) valcritclus
				valcritclus=valcritclus/100D0
			else if (isel2==2) then
				write(*,"(a)") " Input a value, e.g. 6"
				write(*,"(a)") " Note: Inputting P means for any two attractors that have relative value difference less than the one set by option 1, &
				if their interval is less than P*sqrt(dx^2+dy^2+dz^2), &
				where dx,dy,dz are grid spacing in X,Y,Z,then they will be clustered together. If you want to nullify the clustering step after generating &
				basins, simply set this value to 0"
				read(*,*) mergeattdist
			else if (isel2==3) then
				if (numatt<2) then
					write(*,*) "Error: At least two attractors must be presented!"
					cycle
				end if
				write(*,*) "Input the attractors you want to cluster, e.g. 4,5,8,9,22,21"
				read(*,"(a)") c2000tmp
				call str2arr(c2000tmp,ntobeclus) !Find how many terms
				if (allocated(attconvold)) deallocate(attconvold,realattconv,usercluslist)
				allocate(attconvold(-2:numatt),realattconv(-2:numrealatt),usercluslist(ntobeclus))
				realattconv(-2)=-2
				realattconv(-1)=-1
				realattconv(0)=0
				call str2arr(c2000tmp,ntobeclus,usercluslist)
				call sort(usercluslist)
				attconvold=attconv
				idesrealatt=usercluslist(1)
				do itmp=2,ntobeclus
					irealatt=usercluslist(itmp)
					do jtmp=1,nrealatthas(irealatt)
						iatt=realatttable(irealatt,jtmp)
						attconv(iatt)=idesrealatt
					end do
				end do
				call clusdegenatt(1)
				realattconv(1)=1
				do itmp=1,numatt !Build conversion relationship between previous real attractors and new real attractors
					if (itmp/=1) then
						if (attconvold(itmp)/=attconvold(itmp-1)) realattconv(attconvold(itmp))=attconv(itmp)
					end if
				end do
				do iz=2,nz-1
					do iy=2,ny-1
						do ix=2,nx-1
							gridbas(ix,iy,iz)=realattconv(gridbas(ix,iy,iz))
						end do
					end do
				end do
				call detectinterbasgrd(6) !Detect interbasin grids
				write(*,*) "Done!"
				write(*,*)
			end if
		end do
	
	else if (isel==-5) then !Output cube file for basins for visualizing interbasin surface
		write(*,"(a)") " Inputting e.g. 4,6-9,15,17 will output these basins to individual basin[index].cub files to current folder"
		write(*,"(a)") " Inputting ""a"" will output all basins into basin.cub in current folder, the grid value corresponds to basin index"
		read(*,"(a)") c2000tmp
		if (index(c2000tmp,'a')==0) then
            call str2arr(c2000tmp,ntmp)
            if (allocated(tmparr)) deallocate(tmparr)
            allocate(tmparr(ntmp))
            call str2arr(c2000tmp,ntmp,tmparr)
            if (any(tmparr<=0).or.any(tmparr>numatt)) then
                write(*,*) "Error: One or more basin indices exceed valid range!"
                write(*,*) "Press ENTER button to continue"
                read(*,*)
                cycle
            end if
            write(*,"(a)") " Select how to output cube file. The grid value in the select and unselected region will be 1.0 and 0.0, respectively"
            write(*,*) "1 Output all basin grids"
            write(*,*) "2 Output boundary basin grids"
            write(*,*) "3 Output all basin grids where electron density > 0.001 a.u."
		    read(*,*) igrid
		    if (allocated(cubmattmp)) deallocate(cubmattmp)
		    allocate(cubmattmp(nx,ny,nz))
		    do idx=1,ntmp
                ibas=tmparr(idx)
			    write(basinfilename,"('basin',i4.4,'.cub')") ibas
			    write(*,"(' Outputting basin',i8,' as ',a)") ibas,trim(basinfilename)
			    cubmattmp=gridbas(:,:,:)
			    where(cubmattmp/=ibas)
				    cubmattmp=0
			    elsewhere
				    cubmattmp=1
			    end where
			    if (igrid==1) then
                    continue
                else if (igrid==2) then
                    where (.not.interbasgrid) cubmattmp=0
                else if (igrid==3) then
                    if (ifuncbasin==1) then !The cubmat already records electron density
                        where(cubmat<0.001D0) cubmattmp=0
                    else
                        call saverhocub
                        where(rhocub<0.001D0) cubmattmp=0
                    end if
                end if
			    open(10,file=basinfilename,status="replace")
                call outcube(cubmattmp,nx,ny,nz,orgx,orgy,orgz,gridvec1,gridvec2,gridvec3,10)
			    close(10)
		    end do
		    deallocate(cubmattmp)
			write(*,"(a)") " Done! Cube files for the selected basins have been outputted to current folder"
            if (selectyn=='y'.or.selectyn=='Y') then
			    write(*,"(a)") " Values 1 and 0 in these files indicate that the corresponding point belongs and not belongs to the basin, & 
			    respectively. You can plot isosurface with isovalue=0.5 to visualize basin region"
            else
			    write(*,"(a)") " You can plot isosurface with isovalue=0.5 to visualize the basin surface"
            end if
		else
		    if (allocated(cubmattmp)) deallocate(cubmattmp)
		    allocate(cubmattmp(nx-2,ny-2,nz-2))
            cubmattmp=gridbas(2:nx-1,2:ny-1,2:nz-1)
            open(10,file="basin.cub",status="replace")
            call outcube(cubmattmp,nx-2,ny-2,nz-2,orgx+dx,orgy+dy,orgz+dz,gridvec1,gridvec2,gridvec3,10)
            close(10)
		    deallocate(cubmattmp)
			write(*,"(a)") " Done! basin.cub has been outputted to current folder. The grid values correspond to basin index"
		end if
		
	else if (isel==-4) then
        write(*,*) "0 Return"
        write(*,*) "1 Export coordinates of all attractors as attractors.pdb"
        write(*,"(a)") " 2 Export coordinates and function values of all attractors as attractors.pqr"
        write(*,"(a)") " 3 Export coordinates and function values of all attractors as attractors.txt"
        read(*,*) isel2
        if (isel2==0) then
            continue
        else if (isel2==1) then
		    open(10,file="attractors.pdb",status="replace")
		    do iatt=1,numatt
			    write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") "HETATM",iatt,' '//"C "//' ',"ATT",'A',attconv(iatt),attxyz(iatt,:)*b2a,1.0,0.0,"C "
		    end do
		    close(10)
		    write(*,*) "Done, all attractors have been exported to attractors.pdb in current folder"
		    write(*,"(a)") " Note: The residue indices in the pdb file correspond to the attractor indices after clustering, &
            while the atomic indices correspond to the raw attractor indices"
        else if (isel2==2) then
		    open(10,file="attractors.pqr",status="replace")
		    do iatt=1,numatt
			    write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,E18.8E3,f5.2,1x,a2)") "HETATM",iatt,' '//"C "//' ',"ATT",'A',attconv(iatt),attxyz(iatt,:)*b2a,attval(iatt),1.0,"C "
		    end do
		    close(10)
		    write(*,*) "Done, all attractors have been exported to attractors.pqr in current folder"
		    write(*,"(a)") " Note: The residue indices in the pqr file correspond to the attractor indices after clustering, &
            while the atomic indices correspond to the raw attractor indices. The atomic charge column corresponds to function value"
        else if (isel2==3) then
		    open(10,file="attractors.txt",status="replace")
		    do iatt=1,numatt
                write(10,"(3f12.6,1PE18.8E3)") attxyz(iatt,:),attval(iatt)
            end do
		    close(10)
		    write(*,*) "Done, all attractors have been exported to attractors.txt in current folder"
            write(*,"(a)") " In this file, the first three columns correspond to X,Y,Z coordinate in Bohr, the final column is function value"
        end if
		
	else if (isel==-3) then
		write(*,*) "  Attractor       X,Y,Z coordinate (Angstrom)            Value"
		do irealatt=1,numrealatt
			write(*,"(i8,3f14.8,1E18.9)") irealatt,realattxyz(irealatt,1:3)*b2a,realattval(irealatt)
		end do
		do irealatt=1,numrealatt
			if (nrealatthas(irealatt)/=1) then
				write(*,*)
				write(*,"(' The members of degenerate attractor',i6,':')") irealatt
				do itmp=1,nrealatthas(irealatt)
					iatt=realatttable(irealatt,itmp)
					write(*,"(i8,'  XYZ:',3f13.7,'  Value:',1E18.9)") iatt,attxyz(iatt,:)*b2a,attval(iatt)
				end do
			end if
		end do
        if (numrealatt>1) then
            write(*,"(/,a)") " Do you also want to print the attractor information after sorting according to values? (y/n)"
            read(*,*) selectyn
            if (selectyn=='y'.or.selectyn=='Y') then
                if (allocated(tmparr)) deallocate(tmparr)
                allocate(tmparr(numrealatt))
                forall(i=1:numrealatt) tmparr(i)=i
                allocate(tmparrf(numrealatt))
                tmparrf=realattval
                call sortr8(tmparrf,"val",tmparr)
		        write(*,*) "  Attractor       X,Y,Z coordinate (Angstrom)            Value"
		        do idx=1,numrealatt
                    irealatt=tmparr(idx)
			        write(*,"(i8,3f14.8,1E18.9)") irealatt,realattxyz(irealatt,1:3)*b2a,realattval(irealatt)
                end do
                deallocate(tmparrf,tmparr)
            end if
        end if
		
	else if (isel==-2) then
		write(*,*) "q = Quit. Selection method:"
		write(*,*) "a? = Atom ?"
		write(*,*) "c? = Attractor ? (If the attractor is degenerate, average coordinate is used)"
		write(*,*) "d? = Attractor ? (Only for degenerate attractors, cycle each of its members)"
		write(*,*) "For example:"
		write(*,*) """a1 c3"" returns the distance between atom1 and att.3"
		write(*,*) """a4 a2"" returns the distance between atom4 and atom2"
		write(*,*) """c6 a2 a5"" returns the angle of att.6-atom2-atom5"
		write(*,*) """c2 c4 a3 c7"" returns the dihedral angle of att.2-att.4-atom3-att.7"
		write(*,*) """d3 a1"" returns the distance between members of att.3 and atom1"
		write(*,"(a)") " ""d3 a1 c2"" returns the the angle of (members of att.3)--atom1--att.2, &
		meanwhile outputs the vertical distance from (members of att.3) to the line linking atom1 and att.2"
		do while(.true.)
			read(*,"(a)") c80tmp
			c80tmp=adjustl(c80tmp)
			imeasure=0
			do ichar=1,len_trim(c80tmp) !imeasure=1/2/3: measure distance,angle,dihedral
				if (c80tmp(ichar:ichar)==','.or.c80tmp(ichar:ichar)==' ') imeasure=imeasure+1
			end do
			nelement=0 !Validate "d" selection
			do ichar=1,len_trim(c80tmp)
				if (c80tmp(ichar:ichar)=='d') nelement=nelement+1
			end do
			if (nelement/=0) then
				if (c80tmp(1:1)/='d'.or.nelement>1) then
					write(*,*) "Error: ""d"" type of selection can only occur once and must be the first term"
					cycle
				else if (imeasure==3) then
					write(*,*) "Error: Dihedral angle calculation does not support ""d"" type of selection"
					cycle					
				end if
			end if
							
			if (c80tmp(1:1)=='q') then
				exit
			else if (imeasure==1.or.imeasure==2.or.imeasure==3) then
				if (imeasure==1) read(c80tmp,*) ctmp1,ctmp2 !Read two terms
				if (imeasure==2) read(c80tmp,*) ctmp1,ctmp2,ctmp3 !Read three terms
				if (imeasure==3) read(c80tmp,*) ctmp1,ctmp2,ctmp3,ctmp4 !Read four terms
				
				if (ctmp1(1:1)=='a') then
					read(ctmp1(2:),*) iatm
					tmpx1=a(iatm)%x
					tmpy1=a(iatm)%y
					tmpz1=a(iatm)%z
				else if (ctmp1(1:1)=='c') then
					read(ctmp1(2:),*) irealatt
					tmpx1=realattxyz(irealatt,1)
					tmpy1=realattxyz(irealatt,2)
					tmpz1=realattxyz(irealatt,3)
				else if (ctmp1(1:1)=='d') then
					read(ctmp1(2:),*) irealattde
				end if
				if (ctmp2(1:1)=='a') then
					read(ctmp2(2:),*) iatm
					tmpx2=a(iatm)%x
					tmpy2=a(iatm)%y
					tmpz2=a(iatm)%z
				else if (ctmp2(1:1)=='c') then
					read(ctmp2(2:),*) irealatt
					tmpx2=realattxyz(irealatt,1)
					tmpy2=realattxyz(irealatt,2)
					tmpz2=realattxyz(irealatt,3)
				end if
				
				if (imeasure==1) then
					if (ctmp1(1:1)/='d') then
						write(*,"(' The distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") &
						dsqrt((tmpx1-tmpx2)**2+(tmpy1-tmpy2)**2+(tmpz1-tmpz2)**2),dsqrt((tmpx1-tmpx2)**2+(tmpy1-tmpy2)**2+(tmpz1-tmpz2)**2)*b2a
					else if (ctmp1(1:1)=='d') then
						distmin=9999999
						distmax=-1
						distavg=0
						do itmp=1,nrealatthas(irealattde)
							iatt=realatttable(irealattde,itmp)
							tmpx1=attxyz(iatt,1)
							tmpy1=attxyz(iatt,2)
							tmpz1=attxyz(iatt,3)
							distnow=dsqrt((tmpx1-tmpx2)**2+(tmpy1-tmpy2)**2+(tmpz1-tmpz2)**2)
							distavg=distavg+distnow
							if (distnow<distmin) distmin=distnow
							if (distnow>distmax) distmax=distnow
						end do
						distavg=distavg/nrealatthas(irealattde)
						write(*,"(' Note: Attractor',i6,' has',i6,' members')") irealattde,nrealatthas(irealattde)
						write(*,"(' The minimum distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") distmin,distmin*b2a
						write(*,"(' The maximum distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") distmax,distmax*b2a
						write(*,"(' The average distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") distavg,distavg*b2a
					end if
				end if

				if (imeasure==2.or.imeasure==3) then !Analyze one more term, then print angle
					if (ctmp3(1:1)=='a') then
						read(ctmp3(2:),*) iatm
						tmpx3=a(iatm)%x
						tmpy3=a(iatm)%y
						tmpz3=a(iatm)%z
					else if (ctmp3(1:1)=='c') then
						read(ctmp3(2:),*) irealatt
						tmpx3=realattxyz(irealatt,1)
						tmpy3=realattxyz(irealatt,2)
						tmpz3=realattxyz(irealatt,3)
					end if
				end if
				if (imeasure==2) then
					if (ctmp1(1:1)/='d') then
						write(*,"(' The angle is',f12.6,' degree')") xyz2angle(tmpx1,tmpy1,tmpz1,tmpx2,tmpy2,tmpz2,tmpx3,tmpy3,tmpz3)
					else if (ctmp1(1:1)=='d') then
						angmin=180
						angmax=-1
						angavg=0
						distmin=999999
						distmax=-1
						distavg=0
						prjx=0
						prjy=0
						prjz=0
						do itmp=1,nrealatthas(irealattde)
							iatt=realatttable(irealattde,itmp)
							tmpx1=attxyz(iatt,1)
							tmpy1=attxyz(iatt,2)
							tmpz1=attxyz(iatt,3)
							angnow=xyz2angle(tmpx1,tmpy1,tmpz1,tmpx2,tmpy2,tmpz2,tmpx3,tmpy3,tmpz3)
							angavg=angavg+angnow
							if (angnow<angmin) angmin=angnow
							if (angnow>angmax) angmax=angnow
							distnow=potlinedis(tmpx1,tmpy1,tmpz1,tmpx2,tmpy2,tmpz2,tmpx3,tmpy3,tmpz3)
							distavg=distavg+distnow
							if (distnow<distmin) distmin=distnow
							if (distnow>distmax) distmax=distnow
							call pointprjline(tmpx1,tmpy1,tmpz1,tmpx2,tmpy2,tmpz2,tmpx3,tmpy3,tmpz3,prjxtmp,prjytmp,prjztmp)
							prjx=prjx+prjxtmp
							prjy=prjy+prjytmp
							prjz=prjz+prjztmp
						end do
						angavg=angavg/nrealatthas(irealattde)
						distavg=distavg/nrealatthas(irealattde)
						prjx=prjx/nrealatthas(irealattde)
						prjy=prjy/nrealatthas(irealattde)
						prjz=prjz/nrealatthas(irealattde)
						write(*,"(' Note: Attractor',i6,' has',i6,' members')") irealattde,nrealatthas(irealattde)
						write(*,"(' The minimum angle is',f12.6,' degree')") angmin
						write(*,"(' The maximum angle is',f12.6,' degree')") angmax
						write(*,"(' The average angle is',f12.6,' degree')") angavg
						write(*,"(' The minimum distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") distmin,distmin*b2a
						write(*,"(' The maximum distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") distmax,distmax*b2a
						write(*,"(' The average distance is',f12.6,' Bohr (',f12.6 ' Angstrom)')") distavg,distavg*b2a
						write(*,"(' The average X,Y,Z coordinate by projecting the members of ',a,' to the line linking ',a,' and ',a)") trim(ctmp1),trim(ctmp2),trim(ctmp3)
						write(*,"(3f14.8,'   Angstrom')") prjx*b2a,prjy*b2a,prjz*b2a
					end if
				end if
				
				if (imeasure==3) then !Analyze one more term, then print dihedral angle
					if (ctmp4(1:1)=='a') then
						read(ctmp4(2:),*) iatm
						tmpx4=a(iatm)%x
						tmpy4=a(iatm)%y
						tmpz4=a(iatm)%z
					else if (ctmp4(1:1)=='c') then
						read(ctmp4(2:),*) irealatt
						tmpx4=realattxyz(irealatt,1)
						tmpy4=realattxyz(irealatt,2)
						tmpz4=realattxyz(irealatt,3)
					end if
					write(*,"(' The dihedral angle is',f12.6,' degree')") xyz2dih(tmpx1,tmpy1,tmpz1,tmpx2,tmpy2,tmpz2,tmpx3,tmpy3,tmpz3,tmpx4,tmpy4,tmpz4)
				end if
			else
				write(*,*) "Error: Invalid input"
			end if
		end do
		
	else if (isel==-1) then
		write(*,*) "1: On-grid method, Comput .Mat. Sci., 36, 354 (2006)"
		write(*,*) "2: Near-grid method, J. Phys.: Condens. Matter, 21, 08420 (2009)"
		write(*,*) "3: Near-grid method with boundary refinement step"
		write(*,"(a)") " Note: Near-grid method (adapted by Tian Lu) is more accurate than On-grid method and thus is more recommended; with the boundary refinement step, the result will be better"
		read(*,*) igridmethod
		
	else if (isel==0) then
		ioldtextheigh=textheigh
		textheigh=40 !Default textheigh is too small
		call drawbasinintgui
		textheigh=ioldtextheigh
		
	else if (isel==1) then
		isourcedata=1 !Default status is calculating grid data here
		if (allocated(cubmat)) then
			write(*,"(a)") " Note: There has been a grid data in the memory, please select generating the basins by which manner"
			write(*,*) "0 Return"
			write(*,*) "1 Generate the basins by selecting a real space function"
			write(*,*) "2 Generate the basins by using the grid data stored in memory"
			read(*,*) isourcedata
		end if
		if (isourcedata==0) then
			cycle
		else if (isourcedata==1) then
			write(*,*) "Select the real space function to be integrated"
			call selfunc_interface(1,ifuncbasin)
			call setgridforbasin(ifuncbasin)
			if (allocated(cubmat)) deallocate(cubmat)
			allocate(cubmat(nx,ny,nz))
			call savecubmat(ifuncbasin,0,iorbsel)
		else if (isourcedata==2) then
			ifuncbasin=1 !I assume that the file loaded records electron density
			write(*,"(' The range of x is from ',f12.6,' to ',f12.6,' Bohr,' i5,' points')") ,orgx,orgx+(nx-1)*dx,nx
			write(*,"(' The range of y is from ',f12.6,' to ',f12.6,' Bohr,',i5,' points')") ,orgy,orgy+(ny-1)*dy,ny
			write(*,"(' The range of z is from ',f12.6,' to ',f12.6,' Bohr,',i5,' points')") ,orgz,orgz+(nz-1)*dz,nz
			write(*,"(' Total number of grid points is ',i10)") nx*ny*nz
			write(*,"(' Grid spacing in X,Y,Z (Bohr):',3f12.6)") dx,dy,dz
		end if
		
		allocate(grdposneg(nx,ny,nz)) !Record which grids have negative value
		grdposneg=.true.
		do iz=1,nz
			do iy=1,ny
				do ix=1,nx
					if (cubmat(ix,iy,iz)<0D0) then
						grdposneg(ix,iy,iz)=.false.
						cubmat(ix,iy,iz)=-cubmat(ix,iy,iz)
					end if
				end do
			end do
		end do
! 		where (cubmat<0D0)  !DO NOT USE THIS, because I found that "where" will consuming vary large amount of memory!
! 			grdposneg=.false.
! 			cubmat=-cubmat !Invert negative values to positive, after basins are generated the values will be recovered
! 		end where
		if (allocated(gridbas)) deallocate(gridbas)
		allocate(gridbas(nx,ny,nz))
		call setupmovevec
		write(*,*)
		write(*,*) "Generating basins, please wait..."
		call walltime(walltime1)
		call generatebasin(igridmethod) !Generate basins
		call walltime(walltime2)
		write(*,"(' Generating basins took up wall clock time',i10,' s')") walltime2-walltime1	
		numunassign=count(gridbas(2:nx-1,2:ny-1,2:nz-1)==0)
		write(*,"(' The number of unassigned grids:',i12)") numunassign
		numgotobound=count(gridbas(2:nx-1,2:ny-1,2:nz-1)==-1)
		write(*,"(' The number of grids travelled to box boundary:',i12)") numgotobound
		where (.not.grdposneg) cubmat=-cubmat !Recover original grid data
		deallocate(grdposneg)
		
		do iatt=1,numatt !Eliminate the attractors with very low value
			if ( abs(cubmat(attgrid(iatt,1),attgrid(iatt,2),attgrid(iatt,3)))<threslowvalatt ) then
				write(*,"(' Note: There are attractors having very low absolute value (<',1PE8.2,') and thus insignificant, how to deal with them?')") threslowvalatt
				write(*,*) "1 Do nothing"
				write(*,*) "2 Set corresponding grids as unassigned status"
				write(*,*) "3 Assign corresponding grids to the nearest significant attractors"
				write(*,*) "Hint: For most cases, option 3 is recommended"
				read(*,*) isel2
				call elimlowvalatt(threslowvalatt,isel2)
				exit
			end if
		end do
		
		!Generate actual coordinate of attractors, which will be used to plot in the local GUI, and in any other external GUIs
		if (allocated(attxyz)) deallocate(attxyz,attval)
		allocate(attxyz(numatt,3),attval(numatt))
		do iatt=1,numatt
			ix=attgrid(iatt,1)
			iy=attgrid(iatt,2)
			iz=attgrid(iatt,3)
			attxyz(iatt,1)=orgx+(ix-1)*dx
			attxyz(iatt,2)=orgy+(iy-1)*dy
			attxyz(iatt,3)=orgz+(iz-1)*dz
			attval(iatt)=cubmat(ix,iy,iz)
		end do
		!Currently attractors have been finally determined, one shouldn't perturb them further more
		
		call clusdegenatt(0) !Cluster degenerate attractors as "real attractors" and calculate average coordinate and value for the real attractors

		call detectinterbasgrd(6) !Detect interbasin grids
		numinterbas=count(interbasgrid)
		write(*,"(' The number of interbasin grids:',i12)") numinterbas
		
	else if (isel==2) then
		call integratebasin
		
	else if (isel==3.or.isel==4.or.isel==5.or.isel==6.or.isel==7.or.isel==-7.or.isel==8) then
		if (.not.allocated(b)) then
			write(*,"(a)") " Note: No GTF (Gaussian type function) information is available in your input file, please input the file &
			containing GTF information of your system, such as .wfn/.wfx and .fch file. e.g. C:\abc.wfn"
			read(*,"(a)") c200tmp
			call readinfile(c200tmp,1)
		end if
		if (isel==3) then
			call multipolebasin
		else if (isel==4) then
			call LIDIbasin(0)
		else if (isel==5) then
			call LIDIbasin(1)
		else if (isel==6) then
			call LIDIbasin(2)
		else if (isel==7) then
			write(*,*) "0 Return"
			write(*,*) "1 Integrate a specific function with atomic-center + uniform grids"
			write(*,*) "2 The same as 1, but with exact refinement of basin boundary"
			write(*,*) "3 The same as 2, but with approximate refinement of basin boundary"
			write(*,*) "Hint:"
			write(*,*) "Accuracy: 2>=3>>1     Time spent: 2>3>>1     Memory requirement: 3>2=1"
	! 		write(*,*) "Robost: 1=2>3"
			read(*,*) iseltmp
			call integratebasinmix(iseltmp)
		else if (isel==-7) then
			call integratebasinmix_LSB
		else if (isel==8) then
			call integratebasinmix(10)
		end if
		
! 	else if (isel==100) then
! 		call integratebasinrefine
	else if (isel==9) then
		call atmpopinbasin
    else if (isel==10) then
        call HELP_HELV
    else if (isel==11) then
        call basinorbcomp
	end if
end do

end subroutine



!!---- setup move vector and move length for the 26 neighbour
subroutine setupmovevec
use defvar
use basinintmod
vec26x=0
vec26y=0
vec26z=0
!The nearest neighbours:
len26(1:2)=dx
vec26x(1)=1
vec26x(2)=-1
len26(3:4)=dy
vec26y(3)=1
vec26y(4)=-1
len26(5:6)=dz
vec26z(5)=1
vec26z(6)=-1
!On the edges:
len26(7:10)=dsqrt(dx*dx+dy*dy)
vec26x(7)=1
vec26y(7)=1
vec26x(8)=-1
vec26y(8)=1
vec26x(9)=-1
vec26y(9)=-1
vec26x(10)=1
vec26y(10)=-1
len26(11:14)=dsqrt(dx*dx+dz*dz)
vec26x(11)=1
vec26z(11)=1
vec26x(12)=-1
vec26z(12)=1
vec26x(13)=-1
vec26z(13)=-1
vec26x(14)=1
vec26z(14)=-1
len26(15:18)=dsqrt(dy*dy+dz*dz)
vec26y(15)=1
vec26z(15)=1
vec26y(16)=1
vec26z(16)=-1
vec26y(17)=-1
vec26z(17)=-1
vec26y(18)=-1
vec26z(18)=1
!At the vertices:
len26(19:26)=dsqrt(dx*dx+dy*dy+dz*dz)
vec26z(19:22)=1
vec26x(19)=1
vec26y(19)=1
vec26x(20)=-1
vec26y(20)=1
vec26x(21)=-1
vec26y(21)=-1
vec26x(22)=1
vec26y(22)=-1
vec26z(23:26)=-1
vec26x(23)=1
vec26y(23)=-1
vec26x(24)=1
vec26y(24)=1
vec26x(25)=-1
vec26y(25)=1
vec26x(26)=-1
vec26y(26)=-1
if (allocated(gridbas)) then
	gridbas=0 !Unassigned state
	gridbas(1,:,:)=-2 !Set status for box boundary grids
	gridbas(nx,:,:)=-2
	gridbas(:,1,:)=-2
	gridbas(:,ny,:)=-2
	gridbas(:,:,1)=-2
	gridbas(:,:,nz)=-2
end if
end subroutine



!!------- Generate basins from regular grid
! igridmethod=1: On grid, Comput.Mat.Sci.,36,354    =2: Near-grid method (slower, but more accurate), see J.Phys.:Condens.Matter,21,084204
! =3: Near-grid with refinement
! The near-grid method is improved by Tian Lu, namely at later stage automatically switch to on-grid method to guarantee convergence
subroutine generatebasin(igridmethod)
use defvar
use basinintmod
implicit real*8(a-h,o-z)
integer,parameter :: nmaxtrjgrid=3000
integer igridmethod
integer ntrjgrid !Recording trjgrid contains how many elements now
integer trjgrid(nmaxtrjgrid,3) !The trajectory contains which grids, record their indices sequentially, trjgrid(i,1/2/3) = ix/iy/iz of the ith grid
! real*8 trjval(nmaxtrjgrid),gradmaxval(nmaxtrjgrid) !******For debugging******
if (allocated(attgrid)) deallocate(attgrid)
allocate(attgrid(nint(nx*ny*nz/20D0),3)) !I think the number of attractors in general is impossible to exceeds nx*ny*nz/20
numatt=0
write(*,*) "  Attractor       X,Y,Z coordinate (Angstrom)                Value"
!Cycle all grids, but the box boundary ones are ignored (due to gridbas=-2), since can't evalute its gradients in all directions by finite difference

nsteplimit=min( nmaxtrjgrid,nint(dsqrt(dfloat(nx*nx+ny*ny+nz*nz))*2) )
nstepdiscorr=nint(nsteplimit/2D0)
if (igridmethod==1) then
1	continue
	!$OMP PARALLEL DO private(ix,iy,iz,ntrjgrid,inowx,inowy,inowz,trjgrid,valnow,imove,gradtmp,igradmax,gradmax,iatt,itrjgrid,idtmp) &
	!$OMP shared(gridbas,numatt,attgrid) schedule(DYNAMIC) NUM_THREADS(nthreads)
	do iz=2,nz-1
		do iy=2,ny-1
			do ix=2,nx-1
				if (gridbas(ix,iy,iz)/=0) cycle
! 				if (interbasgrid(ix,iy,iz)==.false.) cycle
				ntrjgrid=0
				inowx=ix
				inowy=iy
				inowz=iz
				do while(.true.) !Steepest ascent
					ntrjgrid=ntrjgrid+1
					if (ntrjgrid>nsteplimit) exit !Unconverged. The gridbas for these unassigned grids will still be 0
					trjgrid(ntrjgrid,1)=inowx
					trjgrid(ntrjgrid,2)=inowy
					trjgrid(ntrjgrid,3)=inowz
					!Test all 26 directions to find out the maximal gradient direction
					valnow=cubmat(inowx,inowy,inowz)
					do imove=1,26
						gradtmp=(cubmat(inowx+vec26x(imove),inowy+vec26y(imove),inowz+vec26z(imove))-valnow)/len26(imove)
						if (imove==1.or.gradtmp>gradmax) then
							igradmax=imove
							gradmax=gradtmp
						end if
					end do
					!Test if this is an attractor, if yes, assign all grid in this trajectory
					if (gradmax<=0) then !Equal sign is important, because when system has symmetry, adjacent grid may be degenerate about mirrow plane, now the ascent should be terminated
						if (valnow==0D0) exit !The region far beyond system, the value may be exactly zero due to cutoff of exponent, these grids shouldn't be regarded as attractors
						!$OMP CRITICAL
cyciatt3:				do iatt=1,numatt
							if (inowx==attgrid(iatt,1).and.inowy==attgrid(iatt,2).and.inowz==attgrid(iatt,3)) then
								do itrjgrid=1,ntrjgrid
									gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=iatt
								end do
								exit cyciatt3
							end if
							do imove=1,26
								if (inowx+vec26x(imove)==attgrid(iatt,1).and.inowy+vec26y(imove)==attgrid(iatt,2).and.inowz+vec26z(imove)==attgrid(iatt,3)) then
									do itrjgrid=1,ntrjgrid
										gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=iatt
									end do
									exit cyciatt3
								end if
							end do
						end do cyciatt3
						if (iatt>numatt) then !A new attractor, iatt=numatt+1 currently
							numatt=numatt+1
							do itrjgrid=1,ntrjgrid
								gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=numatt
							end do
							attgrid(numatt,1)=inowx
							attgrid(numatt,2)=inowy
							attgrid(numatt,3)=inowz
							if (grdposneg(inowx,inowy,inowz)) then
								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,cubmat(inowx,inowy,inowz)
							else !This grid should has negative value
								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,-cubmat(inowx,inowy,inowz)
							end if
						end if
						!$OMP end CRITICAL
						exit
					end if
					!Move to next grid
					inowx=inowx+vec26x(igradmax)
					inowy=inowy+vec26y(igradmax)
					inowz=inowz+vec26z(igradmax)
					!Test if this grid has already been assigned
					idtmp=gridbas(inowx,inowy,inowz)
					if (idtmp>0) then
						do itrjgrid=1,ntrjgrid
							gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=idtmp
						end do
						exit
					end if
					!Test if encountered box boundary
					if (inowx==1.or.inowx==nx.or.inowy==1.or.inowy==ny.or.inowz==1.or.inowz==nz) then
						do itrjgrid=1,ntrjgrid
							gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=-1
						end do
						exit
					end if
				end do !End ascent
			end do !End cycle x
		end do !End cycle y
	end do !End cycle z
	!$OMP END PARALLEL DO
	
else if (igridmethod==2.or.igridmethod==3) then
	!$OMP PARALLEL DO private(ix,iy,iz,corrx,corry,corrz,ntrjgrid,inowx,inowy,inowz,trjgrid,valnow,imove,gradtmp,igradmax,gradmax,iatt,&
	!$OMP itrjgrid,idtmp,icorrx,icorry,icorrz,gradx,grady,gradz,sclgrad,ineiidx) shared(gridbas,numatt,attgrid) schedule(DYNAMIC) NUM_THREADS(nthreads)
	do iz=2,nz-1
		do iy=2,ny-1
			do ix=2,nx-1
				if (gridbas(ix,iy,iz)/=0) cycle
				ntrjgrid=0
				inowx=ix
				inowy=iy
				inowz=iz
				corrx=0D0 !Correction vector
				corry=0D0
				corrz=0D0
				do while(.true.) !Steepest ascent
					ntrjgrid=ntrjgrid+1
					trjgrid(ntrjgrid,1)=inowx
					trjgrid(ntrjgrid,2)=inowy
					trjgrid(ntrjgrid,3)=inowz
					if (ntrjgrid>nsteplimit) then !These gridbas for these unconverged grids will still be 0
! 						do itmp=1,ntrjgrid-1 !******For debugging******
! 							write(*,"(i5,3i6,2f16.10)") itmp,trjgrid(itmp,1:3),trjval(itmp),gradmaxval(itmp)
! 						end do
						exit
					end if
					!Test all 26 directions to find out the maximal gradient direction
					valnow=cubmat(inowx,inowy,inowz)
					do imove=1,26
						gradtmp=(cubmat(inowx+vec26x(imove),inowy+vec26y(imove),inowz+vec26z(imove))-valnow)/len26(imove)
						if (imove==1.or.gradtmp>gradmax) then
							igradmax=imove
							gradmax=gradtmp
						end if
					end do
!  					write(*,"(3i5,2f14.8)") inowx,inowy,inowz,gradmax,valnow !Trace the trajectory !******For debugging******
					if (gradmax<=0D0) then !Equal sign is important, because when system has symmetry, adjacent grid may be degenerate about mirrow plane, now the ascent should be terminated
						if (valnow==0D0) exit !The region far beyond system, the value may be exactly zero due to cutoff of exponent, these grids shouldn't be regarded as attractors
						!$OMP CRITICAL
cyciatt:				do iatt=1,numatt
							!Test if current grid is attractor iatt, is yes, assign all grid in this trajectory
							if (inowx==attgrid(iatt,1).and.inowy==attgrid(iatt,2).and.inowz==attgrid(iatt,3)) then
								do itrjgrid=1,ntrjgrid
									gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=iatt
								end do
								exit cyciatt
							end if
							!Test if neighbour grid (+/-x,+/-y,+/-z) is attractor iatt, is yes, assign all grid in this trajectory
							!The reason I do this is because when the points are symmetric to Cartesian plane, many adjacent and degenerate attractors will occur,&
							!while this treatment can combine them as a single one
							do imove=1,26
								if (inowx+vec26x(imove)==attgrid(iatt,1).and.inowy+vec26y(imove)==attgrid(iatt,2).and.inowz+vec26z(imove)==attgrid(iatt,3)) then
									do itrjgrid=1,ntrjgrid
										gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=iatt
									end do
									exit cyciatt
								end if
							end do
						end do cyciatt
						!A new attractor
						if (iatt==numatt+1) then
							numatt=numatt+1
							do itrjgrid=1,ntrjgrid
								gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=numatt
							end do
							attgrid(numatt,1)=inowx
							attgrid(numatt,2)=inowy
							attgrid(numatt,3)=inowz
							if (grdposneg(inowx,inowy,inowz)) then
								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,cubmat(inowx,inowy,inowz)
							else !This grid should has negative value
								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,-cubmat(inowx,inowy,inowz)
							end if
						end if
						!$OMP end CRITICAL
						exit
					end if
					!Correction step may lead to oscillator when encountering circular ELF/LOL attractor, so if the current number of step is already large
					!(larger than half of upper limit of step number), then correction will be disabled (namely switch to on-grid method) to guarantee convergence.
					if ( ntrjgrid<nstepdiscorr .and. (abs(corrx)>(dx/2D0).or.abs(corry)>(dy/2D0).or.abs(corrz)>(dz/2D0)) ) then !This time we do correction step
						if (abs(corrx)>(dx/2D0)) then
							icorrx=nint(corrx/abs(corrx)) !Get sign of corrx
							inowx=inowx+icorrx
							corrx=corrx-icorrx*dx
						end if
						if (abs(corry)>(dy/2D0)) then
							icorry=nint(corry/abs(corry))
							inowy=inowy+icorry
							corry=corry-icorry*dy
						end if
						if (abs(corrz)>(dz/2D0)) then
							icorrz=nint(corrz/abs(corrz))
							inowz=inowz+icorrz
							corrz=corrz-icorrz*dz
						end if
					else !Move to next grid according to maximal gradient and then update correction vector
						!Calculate true gradient
						gradx=(cubmat(inowx+1,inowy,inowz)-cubmat(inowx-1,inowy,inowz))/(2*dx)
						grady=(cubmat(inowx,inowy+1,inowz)-cubmat(inowx,inowy-1,inowz))/(2*dy)
						gradz=(cubmat(inowx,inowy,inowz+1)-cubmat(inowx,inowy,inowz-1))/(2*dz)
						sclgrad=min(dx/abs(gradx),dy/abs(grady),dz/abs(gradz))
						inowx=inowx+vec26x(igradmax)
						inowy=inowy+vec26y(igradmax)
						inowz=inowz+vec26z(igradmax)
						corrx=corrx+gradx*sclgrad-vec26x(igradmax)*dx
						corry=corry+grady*sclgrad-vec26y(igradmax)*dy
						corrz=corrz+gradz*sclgrad-vec26z(igradmax)*dz
! 						if (abs(corrx)>(dx/2D0).or.abs(corry)>(dy/2D0).or.abs(corrz)>(dz/2D0)) cycle
					end if
					!Test if this grid has already been assigned and all of the neighbours were also assigned to the same attractor
					idtmp=gridbas(inowx,inowy,inowz)
					if (idtmp>0) then
						do imove=1,26
							ineiidx=gridbas(inowx+vec26x(imove),inowy+vec26y(imove),inowz+vec26z(imove))
							if (ineiidx/=idtmp) exit
						end do
						if (imove==27) then
							do itrjgrid=1,ntrjgrid
								gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=idtmp
							end do
							exit
						end if
					end if
					!Test if encountered box boundary				
					if (inowx==1.or.inowx==nx.or.inowy==1.or.inowy==ny.or.inowz==1.or.inowz==nz) then
						do itrjgrid=1,ntrjgrid
							gridbas(trjgrid(itrjgrid,1),trjgrid(itrjgrid,2),trjgrid(itrjgrid,3))=-1
						end do
						exit
					end if
				end do !End ascent
			end do !End cycle x
		end do !End cycle y
	end do !End cycle z
	!$OMP END PARALLEL DO
	
	!!!!   Refining the basin boundary
	if (igridmethod==3) then
! 		do itime=1,3 !Refine one time in general is sufficient
		write(*,*) "Detecting boundary grids..."
		call detectinterbasgrd(6) !It seems that using 26 directions to determine boundary grids doesn't bring evident benefit
		write(*,"(' There are',i12,' grids at basin boundary')") count(interbasgrid)
		write(*,*) "Refining basin boundary..."
		!Below code is the adapted copy of above near-grid code
		!$OMP PARALLEL DO private(ix,iy,iz,corrx,corry,corrz,ntrjgrid,inowx,inowy,inowz,valnow,imove,gradtmp,igradmax,gradmax,iatt,&
		!$OMP ineiidx,idtmp,icorrx,icorry,icorrz,gradx,grady,gradz,sclgrad) shared(gridbas,attgrid) schedule(DYNAMIC) NUM_THREADS(nthreads)
		do iz=2,nz-1
			do iy=2,ny-1
				do ix=2,nx-1
					if (.not.interbasgrid(ix,iy,iz)) cycle
					if (gridbas(ix,iy,iz)<=0) cycle !Ignored the ones unassigned or gone to box boundary
					ntrjgrid=0
					inowx=ix
					inowy=iy
					inowz=iz
					corrx=0D0 !Correction vector
					corry=0D0
					corrz=0D0
					do while(.true.) !Steepest ascent
						ntrjgrid=ntrjgrid+1
						if (ntrjgrid>nsteplimit) exit
						valnow=cubmat(inowx,inowy,inowz)
						do imove=1,26
							gradtmp=(cubmat(inowx+vec26x(imove),inowy+vec26y(imove),inowz+vec26z(imove))-valnow)/len26(imove)
							if (imove==1.or.gradtmp>gradmax) then
								igradmax=imove
								gradmax=gradtmp
							end if
						end do
						if (gradmax<=0) then !Equal sign is important, because when system has symmetry, adjacent grid may be degenerate about mirrow plane, now the ascent should be terminated
cyciatt2:					do iatt=1,numatt
								if (inowx==attgrid(iatt,1).and.inowy==attgrid(iatt,2).and.inowz==attgrid(iatt,3)) then
									gridbas(ix,iy,iz)=iatt
									exit cyciatt2
								end if
								do imove=1,26 !Test if neighbour grid (+/-x,+/-y,+/-z) is attractor iatt
									if (inowx+vec26x(imove)==attgrid(iatt,1).and.inowy+vec26y(imove)==attgrid(iatt,2).and.inowz+vec26z(imove)==attgrid(iatt,3)) then
										gridbas(ix,iy,iz)=iatt
										exit cyciatt2
									end if
								end do
							end do cyciatt2
							if (iatt>numatt) then
								write(*,*) "Warning: Found new attractor at refining process!"
							end if
							exit
						end if
						if ( ntrjgrid<nstepdiscorr .and. (abs(corrx)>(dx/2D0).or.abs(corry)>(dy/2D0).or.abs(corrz)>(dz/2D0)) ) then !This time we do correction step
							if (abs(corrx)>(dx/2D0)) then
								icorrx=nint(corrx/abs(corrx)) !Get sign of corrx
								inowx=inowx+icorrx
								corrx=corrx-icorrx*dx
							end if
							if (abs(corry)>(dy/2D0)) then
								icorry=nint(corry/abs(corry))
								inowy=inowy+icorry
								corry=corry-icorry*dy
							end if
							if (abs(corrz)>(dz/2D0)) then
								icorrz=nint(corrz/abs(corrz))
								inowz=inowz+icorrz
								corrz=corrz-icorrz*dz
							end if
						else !Move to next grid according to maximal gradient and then update correction vector
							gradx=(cubmat(inowx+1,inowy,inowz)-cubmat(inowx-1,inowy,inowz))/(2*dx)
							grady=(cubmat(inowx,inowy+1,inowz)-cubmat(inowx,inowy-1,inowz))/(2*dy)
							gradz=(cubmat(inowx,inowy,inowz+1)-cubmat(inowx,inowy,inowz-1))/(2*dz)
							sclgrad=min(dx/abs(gradx),dy/abs(grady),dz/abs(gradz))
							inowx=inowx+vec26x(igradmax)
							inowy=inowy+vec26y(igradmax)
							inowz=inowz+vec26z(igradmax)
							corrx=corrx+gradx*sclgrad-vec26x(igradmax)*dx
							corry=corry+grady*sclgrad-vec26y(igradmax)*dy
							corrz=corrz+gradz*sclgrad-vec26z(igradmax)*dz
						end if
						idtmp=gridbas(inowx,inowy,inowz)
						if (ntrjgrid>60.and.idtmp>0) then !If enable this doesn't affect result detectably, but enabling it will evidently reduce computational cost
							do imove=1,26
								ineiidx=gridbas(inowx+vec26x(imove),inowy+vec26y(imove),inowz+vec26z(imove))
								if (ineiidx/=idtmp) exit
							end do
							if (imove==27) then
								gridbas(ix,iy,iz)=idtmp
								exit
							end if
						end if
						if (inowx==1.or.inowx==nx.or.inowy==1.or.inowy==ny.or.inowz==1.or.inowz==nz) exit
					end do !End ascent
				end do !End cycle x
			end do !End cycle y
		end do !End cycle z
		!$OMP END PARALLEL DO
! 		end do
	end if
	
end if
end subroutine



!!------- Eliminate low value attractors and the corresponding basin
! imode=1 Do nothing
! imode=2 Set corresponding grids as unassigned status
! imode=3 Assign corresponding grids to the nearest significant attractors
subroutine elimlowvalatt(threslowvalatt,imode)
use defvar
use basinintmod
implicit real*8(a-h,o-z)
integer imode
real*8 threslowvalatt
integer highvalatt(numatt)
integer att2att(-2:numatt) !Convert indices of old attractors to new ones
integer atttypelist(numatt)
if (imode==1) then
	return
else
	icounthigh=0
	atttypelist=0
	do iatt=1,numatt !Generate list
		if (abs(cubmat(attgrid(iatt,1),attgrid(iatt,2),attgrid(iatt,3)))>=threslowvalatt) then
			icounthigh=icounthigh+1
			highvalatt(icounthigh)=iatt
			atttypelist(iatt)=1 !high
		end if
	end do
	nlowvalatt=numatt-icounthigh
	if (imode==2) then
		do iz=2,nz-1
			do iy=2,ny-1
				do ix=2,nx-1
					if (atttypelist(gridbas(ix,iy,iz))==0) gridbas(ix,iy,iz)=0
				end do
			end do
		end do
	else if (imode==3) then
		do iz=2,nz-1
			rnowz=orgz+(iz-1)*dz
			do iy=2,ny-1
				rnowy=orgy+(iy-1)*dy
				do ix=2,nx-1
					rnowx=orgx+(ix-1)*dx
					if (atttypelist(gridbas(ix,iy,iz))==0) then !Find out which grid belongs to insignificant attractors
						shortmaxsqr=1D20
						do iatt=1,numatt !Will be attribute to the nearest significant attractors
							if (atttypelist(iatt)==1) then
								xatt=orgx+(attgrid(iatt,1)-1)*dx
								yatt=orgy+(attgrid(iatt,2)-1)*dy
								zatt=orgz+(attgrid(iatt,3)-1)*dz
								disttmpsqr=(xatt-rnowx)**2+(yatt-rnowy)**2+(zatt-rnowz)**2
								if (disttmpsqr<shortmaxsqr) then
									ishortmax=iatt
									shortmaxsqr=disttmpsqr
								end if
							end if
						end do
						gridbas(ix,iy,iz)=ishortmax
					end if
				end do
			end do
		end do
	end if
	!Build conversion table between old and new attractors, so that the indices of attractors could be contiguous
	att2att=-3
	att2att(-2)=-2
	att2att(-1)=-1
	att2att(0)=0
	icount=0
	do iatt=1,numatt
		if (atttypelist(iatt)==0) cycle
		icount=icount+1
		att2att(iatt)=icount
	end do
	numatt=numatt-nlowvalatt
	!Update attgrid
	do iatt=1,numatt
		attgrid(iatt,:)=attgrid(highvalatt(iatt),:)
	end do
	!Final update of grid attribution
	do iz=2,nz-1
		do iy=2,ny-1
			do ix=2,nx-1
				gridbas(ix,iy,iz)=att2att(gridbas(ix,iy,iz))
			end do
		end do
	end do
	write(*,"(i6,' insignificant attractors have been eliminated')") nlowvalatt
end if
end subroutine



!!------- !Detect interbasin grids for real attractors
!ndir can be 6 and 26, meaning if only primary 6 directions or all 26 directions are used to determine boundary grids
subroutine detectinterbasgrd(ndir)
use defvar
use basinintmod
implicit real*8(a-h,o-z)
integer ndir
if (allocated(interbasgrid)) deallocate(interbasgrid)
allocate(interbasgrid(nx,ny,nz))
interbasgrid=.false.
do iz=2,nz-1
	do iy=2,ny-1
		do ix=2,nx-1
			idxnow=gridbas(ix,iy,iz)
			do imove=1,ndir !To reduce the number of interbasin grids and thus speed up displaying, I don't use full 26 direction as criterion
				idxmove=gridbas(ix+vec26x(imove),iy+vec26y(imove),iz+vec26z(imove))
				if ((idxmove/=idxnow).and.idxmove/=-2) then !This is a grid between two or more basins, and is not adjacent to box boundary
					interbasgrid(ix,iy,iz)=.true.
					exit
				end if
			end do
		end do
	end do
end do
end subroutine



!!------- Cluster degenerate and adjacent attractors together
!!Nearest neighbour method is used to cluster the attractors, see Leach book p493. Each cluster corresponds to a final attractor
!If imode=0, run the whole subroutine; if imode=1, then external attconv is provided, and only some code in this routine is used to make the index contiguous
subroutine clusdegenatt(imode)
use defvar
use basinintmod
implicit real*8(a-h,o-z)
integer imode
integer neleatt,eleatt(numatt) !eleatt(i)=j means the ith element in clustering stage corresponds to actual jth attractor. nclustlist is its total current elements
integer passlist(numatt) !if the ith element is 1, means i attractor has passed clustering stage
integer ncluster,ncluele(numatt),cluele(numatt,numatt) !cluele(i,4/5/9...) means cluster i has element 4,5,9..., ncluele(i) is current size of cluster i, ncluster is total number of cluster
real*8 attdismat(numatt,numatt),attvalmat(numatt,numatt) !Distance matrix and relative value difference matrix of the original attractors
distcrit=mergeattdist*dsqrt(dx*dx+dy*dy+dz*dz) !Distance criterion
irefined=0 !If 1, that means this combining process takes effects
idebugclusdegenatt=0 !If output debug information

if (imode==1) goto 2

if (allocated(attconv)) deallocate(attconv)
allocate(attconv(-2:numatt))
do iatt=-2,numatt !First assume that each attractor is a real attractor
	attconv(iatt)=iatt
end do

!Generate difference matrix for distance and relative value 
attdismat=0D0
attvalmat=0D0
do iatt=1,numatt
	valiatt=attval(iatt)
	xiatt=attxyz(iatt,1)
	yiatt=attxyz(iatt,2)
	ziatt=attxyz(iatt,3)
	do jatt=iatt+1,numatt
		valjatt=attval(jatt)
		xjatt=attxyz(jatt,1)
		yjatt=attxyz(jatt,2)
		zjatt=attxyz(jatt,3)
		distdiff=dsqrt( (xiatt-xjatt)**2+(yiatt-yjatt)**2+(ziatt-zjatt)**2 )
		attdismat(iatt,jatt)=distdiff
		valdiff=abs((valiatt-valjatt)/valiatt)
		attvalmat(iatt,jatt)=valdiff
		if (valdiff<valcritclus.and.distdiff<distcrit) irefined=1
! 		write(10,"(2i5,f16.10)") iatt,jatt,attvalmat(iatt,jatt)
	end do
end do
attdismat=attdismat+transpose(attdismat)
attvalmat=attvalmat+transpose(attvalmat)
do iatt=1,numatt
	attdismat(iatt,iatt)=attdismat(iatt,iatt)/2D0
	attvalmat(iatt,iatt)=attvalmat(iatt,iatt)/2D0
end do
if (irefined==0) then !Nothing to do. Construct real attractor information and then directly return
	numrealatt=numatt
	if (allocated(nrealatthas)) deallocate(nrealatthas,realattval,realattxyz,realatttable)
	allocate(nrealatthas(numrealatt),realattval(numrealatt),realattxyz(numrealatt,3),realatttable(numrealatt,1))
	nrealatthas=1
	do iatt=1,numatt
		realattxyz(iatt,:)=attxyz(iatt,:)
		realattval(iatt)=attval(iatt)
		realatttable(iatt,1)=iatt
	end do
	return 
else
	write(*,*) "Degenerate attractors detected, clustering them..."
	write(*,"(' Criterion for clustering: Interval <',f8.5,' Bohr, value difference <',f8.5,'%')") distcrit,valcritclus*100
end if
! do i=1,numatt !Print distance matrix
! 	do j=1,numatt
! 		write(10,*) i,j,attdismat(i,j)
! 	end do
! end do

!Clustering via nearest neighbour method, and accordingly update conversion list
passlist=0
do iatt=1,numatt
	if (passlist(iatt)==1) cycle
	if (count(attdismat(iatt,:)<distcrit)==1) then !This attractor has no neighbour
		if (passlist(iatt)==1) cycle
		cycle
	end if
	!Construct mapping table of attractor, capacity is neleatt
	neleatt=0 !The number of degenerate attractors in this batch
	do jatt=iatt,numatt
		if (attvalmat(iatt,jatt)<valcritclus) then
			neleatt=neleatt+1
			eleatt(neleatt)=jatt
			passlist(jatt)=1
		end if
	end do
	if (neleatt==1) cycle !Although this attractor has one or more neighbours, but no other attractor is degenerate with itself, so pass
	if (idebugclusdegenatt==1) then
		write(*,"(a,i5)") " Attractors in this time of clustering triggered by attractor:",iatt
		write(*,"(' The number of initial clusters:',i5,', corresponding to these attractor:')") neleatt
		write(*,"(15i5)") eleatt(1:neleatt)
	end if

	!Start clustering. Now you should forget actual index of attractors. Only consider element 1,2,3...
	!Initialize, assume that each element in this batch of degenerate attractors is a cluster
	ncluster=neleatt
	ncluele(:)=1
	do iclu=1,ncluster
		cluele(iclu,1)=iclu
	end do
	ntimeclu=0
	do while(.true.) !Gradually clustering until closet distance between any cluster pair is larger than criteria
		ntimeclu=ntimeclu+1
		imergeclu=0
		do iclu=1,ncluster !Note that ncluster remain unchanged, eliminated cluster is simply emptyed but still reserve its index
			if (ncluele(iclu)==0) cycle !Has been incorporated to others
			shortmax=9999999D0
			ishortmax=0
			do jclu=iclu+1,ncluster
				if (ncluele(jclu)==0) cycle
				!Calculate shortest distance between cluster i and j
				shortmaxtmp=9999999D0
				do ieletmp=1,ncluele(iclu)
					do jeletmp=1,ncluele(jclu)
						iele=cluele(iclu,ieletmp)
						jele=cluele(jclu,jeletmp)
						distele=attdismat(eleatt(iele),eleatt(jele))
						if (distele<shortmaxtmp) shortmaxtmp=distele
					end do
				end do
				if (shortmaxtmp<shortmax) then !Find out the closest cluster to cluster iclu
					shortmax=shortmaxtmp
					ishortmax=jclu
				end if
			end do
! 			write(*,"(i3,f16.6,4i5,i2)") ntimeclu,shortmax,iclu,ishortmax,abs(shortmax<distcrit)
			if (shortmax<distcrit) then !Combine cluster ishortmax to cluster iclu
				imergeclu=1
				newlen=ncluele(iclu)+ncluele(ishortmax)
				cluele(iclu,ncluele(iclu)+1:newlen)=cluele(ishortmax,1:ncluele(ishortmax))
				ncluele(iclu)=newlen
				ncluele(ishortmax)=0 !Empty
			end if
		end do
		if (imergeclu==0) exit !Between all cluster pair the shortmax is longer than distcrit, therefore exit
	end do
	do iclu=1,ncluster !Update attractor conversion list
		if (ncluele(iclu)==0) cycle
		if (idebugclusdegenatt==1) then
			write(*,"(' Clustering finished after',i5,' times of cycle')") ntimeclu
			write(*,"(' Cluster',i5,', deriving from  attractor',i5,', contains:')") iclu,eleatt(iclu)
	 		write(*,"(' Element:  ',15i5)") cluele(iclu,1:ncluele(iclu))
			write(*,"(' Attractor:',15i5)") eleatt(cluele(iclu,1:ncluele(iclu)))
		end if
		do itmp=1,ncluele(iclu)
			ielimele=cluele(iclu,itmp)
			ielimatt=eleatt(ielimele)
			itargetatt=eleatt(iclu)
			attconv(ielimatt)=itargetatt
		end do
	end do
end do

! do iatt=1,numatt !Print conversion list just after clustering
! 	write(*,*) iatt,attconv(iatt)
! end do
!Update attractor conversion list to make the attractor index contiguous, and hence get "real" attractor
!e.g. old conversion list, after slash is the new one
!   1           1/1
!   2           2/2
!   ...         2/2
!  13           2/2
!  14          14/3 <---After above clustering, 14,16,17 will convert to 14. The first new entry at right column always matched corresponding left column, this is guaranteed by above clustering code
!  15           2/2
!  16          14/3
!  17          14/3
!  18          18/4
2 numrealatt=0
passlist=0 !Record which attractor has already been converted
do iatt=1,numatt
	if (passlist(iatt)==1) cycle
	numrealatt=numrealatt+1 !Of course, numrealatt always <=iatt, so will not unexpectly overlap data during update conversion list
	idestmp=attconv(iatt)
	do jatt=iatt,numatt
		if (attconv(jatt)==idestmp) then
			passlist(jatt)=1
			attconv(jatt)=numrealatt
		end if
	end do
end do
! do iatt=1,numatt !Print conversion list after reordered
! 	write(*,*) iatt,attconv(iatt)
! end do

call updaterealattprop

if (imode==1) return
!Update basin attribution
do iz=2,nz-1
	do iy=2,ny-1
		do ix=2,nx-1
			gridbas(ix,iy,iz)=attconv(gridbas(ix,iy,iz))
		end do
	end do
end do
!Output final attractors
write(*,*) "The attractors after clustering:"
write(*,*) "   Index      Average X,Y,Z coordinate (Angstrom)               Value"
do irealatt=1,numrealatt
	write(*,"(i8,3f15.8,f20.8)") irealatt,realattxyz(irealatt,:)*b2a,realattval(irealatt)
end do
end subroutine



!!------- Update real attractors' properties (average value, xyz, the number of members)
subroutine updaterealattprop
use defvar
use basinintmod
implicit real*8(a-h,o-z)
if (allocated(nrealatthas)) deallocate(nrealatthas,realattval,realattxyz,realatttable)
allocate(nrealatthas(numrealatt),realattval(numrealatt),realattxyz(numrealatt,3),realatttable(numrealatt,numatt))
nrealatthas=0
realattval=0D0
realattxyz=0D0
do iatt=1,numatt
	irealatt=attconv(iatt)
	nrealatthas(irealatt)=nrealatthas(irealatt)+1
	realatttable(irealatt,nrealatthas(irealatt))=iatt
	realattxyz(irealatt,:)=realattxyz(irealatt,:)+attxyz(iatt,:)
	realattval(irealatt)=realattval(irealatt)+attval(iatt)
end do
do irealatt=1,numrealatt
	realattval(irealatt)=realattval(irealatt)/nrealatthas(irealatt)
	realattxyz(irealatt,:)=realattxyz(irealatt,:)/nrealatthas(irealatt)
end do
end subroutine



!!------------------ Set grid for generating basin, adapted from the subroutine "setgrid"
subroutine setgridforbasin(ifuncsel)
use defvar
use GUI
implicit real*8 (a-h,o-z)
integer :: iselexttype=3,ifuncsel
real*8 :: molxlen,molylen,molzlen,tmpx,tmpy,tmpz,rhocrit=1D-6
real*8 :: gridextdist=5D0,enlarbox=2.1D0,spclowqual=0.2D0,spcmedqual=0.1D0,spchighqual=0.06D0,spclunaqual=0.04D0,tmparr6(6)
character c80tmp*80,cubefilename*200
if (.not.allocated(b)) iselexttype=2 !Unable to set box using rho
do while(.true.)
	if (iselexttype==1) then
		orgx=minval(a%x)-gridextdist
		orgy=minval(a%y)-gridextdist
		orgz=minval(a%z)-gridextdist
		endx=maxval(a%x)+gridextdist
		endy=maxval(a%y)+gridextdist
		endz=maxval(a%z)+gridextdist
	else if (iselexttype==2) then
		orgx=minval( a(:)%x-enlarbox*vdwr_tianlu(a(:)%index) )
		orgy=minval( a(:)%y-enlarbox*vdwr_tianlu(a(:)%index) )
		orgz=minval( a(:)%z-enlarbox*vdwr_tianlu(a(:)%index) )
		endx=maxval( a(:)%x+enlarbox*vdwr_tianlu(a(:)%index) )
		endy=maxval( a(:)%y+enlarbox*vdwr_tianlu(a(:)%index) )
		endz=maxval( a(:)%z+enlarbox*vdwr_tianlu(a(:)%index) )
	end if
	molxlen=endx-orgx
	molylen=endy-orgy
	molzlen=endz-orgz
	ntotlow=(nint(molxlen/spclowqual)+1)*(nint(molylen/spclowqual)+1)*(nint(molzlen/spclowqual)+1)
	ntotmed=(nint(molxlen/spcmedqual)+1)*(nint(molylen/spcmedqual)+1)*(nint(molzlen/spcmedqual)+1)
	ntothigh=(nint(molxlen/spchighqual)+1)*(nint(molylen/spchighqual)+1)*(nint(molzlen/spchighqual)+1)
	ntotluna=(nint(molxlen/spclunaqual)+1)*(nint(molylen/spclunaqual)+1)*(nint(molzlen/spclunaqual)+1)
	
	write(*,*) "Please select a method for setting up grid"
	if (iselexttype==1) write(*,"(a,f10.5,a)") " -10 Set grid extension distance for mode 1~6, current: Fixed,",gridextdist," Bohr"
	if (iselexttype==2) write(*,"(a)") " -10 Set grid extension distance for mode 1~6, current: Adaptive"
	if (iselexttype==3) write(*,"(a)") " -10 Set grid extension distance for mode 1~6, current: Detect rho isosurface"
	if (iselexttype==1.or.iselexttype==2) then
		write(*,"(a,f4.2,a,i14)") " 1 Low quality grid, spacing=",spclowqual," Bohr, number of grids:    ",ntotlow
		write(*,"(a,f4.2,a,i14)") " 2 Medium quality grid, spacing=",spcmedqual," Bohr, number of grids: ",ntotmed
		write(*,"(a,f4.2,a,i14)") " 3 High quality grid, spacing=",spchighqual," Bohr, number of grids:   ",ntothigh
		write(*,"(a,f4.2,a,i14)") " 4 Lunatic quality grid, spacing=",spclunaqual," Bohr, number of grids:",ntotluna
	else
		write(*,"(a,f4.2,a,i14)") " 1 Low quality grid, spacing=",spclowqual," Bohr, cost: 1x"
		write(*,"(a,f4.2,a,i14)") " 2 Medium quality grid, spacing=",spcmedqual," Bohr, cost: 8x"
		write(*,"(a,f4.2,a,i14)") " 3 High quality grid, spacing=",spchighqual," Bohr, cost: 36x"
		write(*,"(a,f4.2,a,i14)") " 4 Lunatic quality grid, spacing=",spclunaqual," Bohr, cost: 120x"
	end if
	write(*,*) "5 Only input grid spacing, automatically set other parameters"
	write(*,*) "6 Only input the number of points in X,Y,Z, automatically set other parameters"
	write(*,*) "7 Input original point, translation vector and the number of points"
	write(*,*) "8 Set center position, grid spacing and box length"
	write(*,*) "9 Use grid setting of another cube file"
	write(*,*) "10 Set box of grid data visually using a GUI window"
	
	read(*,*) igridsel
	if (igridsel/=-10) then
		exit
	else
		write(*,*) "Please select the type of extension distance:"
		write(*,*) "1 Use fixed value in each direction"
		write(*,*) "2 Adaptively determined according to van der Waals radius"
		write(*,"(a,1PE7.1)") " 3 Detect rho to make the box just accommodate its isosurface, isoval: ",rhocrit
		write(*,*) "Hint: In general mode 3 is recommended, mode 1 should be avoided to be used"
		read(*,*) iselexttype
		if (iselexttype==1) then
			do while(.true.)
				write(*,*) "Input extension distance (Bohr) e.g. 6.5"
				read(*,*) gridextdist
				if (gridextdist>0D0) exit
				write(*,*) "Error: The value must be larger than 0!" 
			end do
		else if (iselexttype==2) then
			do while(.true.)
				write(*,*) "Input the ratio for scaling vdW radii, e.g. 1.7"
				write(*,"(' Current value is',f12.6)") enlarbox
				read(*,*) enlarbox
				if (enlarbox>0D0) exit
				write(*,*) "Error: The value must be larger than 0!" 
			end do
		else if (iselexttype==3) then
			write(*,*) "Select the isovalue of rho"
			write(*,*) "1: 0.0001       (10^-4)"
			write(*,*) "2: 0.00001      (10^-5)"
			write(*,*) "3: 0.000005   (5*10^-6)"
			write(*,*) "4: 0.000003   (3*10^-6)"
			write(*,*) "5: 0.000001     (10^-6) In general recommended"
			write(*,*) "6: 0.0000005  (5*10^-7)"
			write(*,*) "7: 0.0000001    (10^-7)"
			write(*,*) "8: 0.00000005 (5*10^-8)"
			write(*,*) "9: 0.00000001   (10^-8)"
			write(*,*) "10: Input by yourself"
			read(*,*) iselrho
			if (iselrho==1) rhocrit=1D-4
			if (iselrho==2) rhocrit=1D-5
			if (iselrho==3) rhocrit=5D-6
			if (iselrho==4) rhocrit=3D-6
			if (iselrho==5) rhocrit=1D-6
			if (iselrho==6) rhocrit=5D-7
			if (iselrho==7) rhocrit=1D-7
			if (iselrho==8) rhocrit=5D-8
			if (iselrho==9) rhocrit=1D-8
			if (iselrho==10) then
				write(*,*) "Input the value, e.g. 0.000005 or 5D-6"
				read(*,*) rhocrit
			end if
		end if
	end if
end do

if (igridsel==1.or.igridsel==2.or.igridsel==3.or.igridsel==4.or.igridsel==5.or.igridsel==6) then
	if (iselexttype==1) then
		aug3D=gridextdist !Set aug3D, which only affects GUI display, so that the scope of the axis is wide enough to contain the whole grid region
	else if (iselexttype==2) then
		aug3D=maxval(enlarbox*vdwr_tianlu(a(:)%index))*1.35D0
	else if (iselexttype==3) then
		call detectboxrho(rhocrit) !determine orgx/y/z and endx/y/z
		molxlen=endx-orgx
		molylen=endy-orgy
		molzlen=endz-orgz
		tmparr6(1)=abs(orgx-minval(a%x))
		tmparr6(2)=abs(orgy-minval(a%y))
		tmparr6(3)=abs(orgz-minval(a%z))
		tmparr6(4)=abs(endx-maxval(a%x))
		tmparr6(5)=abs(endy-maxval(a%y))
		tmparr6(6)=abs(endz-maxval(a%z))
		aug3D=maxval(tmparr6)+1
	end if
end if

!Note: orgx,orgy,orgz,endx,endy,endz as well as molx/y/zlen for igridsel==1~6 have already been set above
if (igridsel==1.or.igridsel==2.or.igridsel==3.or.igridsel==4.or.igridsel==5) then
	if (igridsel==1) dx=spclowqual
	if (igridsel==2) dx=spcmedqual
	if (igridsel==3) dx=spchighqual
	if (igridsel==4) dx=spclunaqual
	if (igridsel==5) then
		write(*,*) "Input the grid spacing (bohr)  e.g. 0.08"
		read(*,*) dx
	end if
	dy=dx
	dz=dx
	nx=nint(molxlen/dx)+1
	ny=nint(molylen/dy)+1
	nz=nint(molzlen/dz)+1
else if (igridsel==6) then
	write(*,*) "Input the number of grid points in X,Y,Z direction   e.g. 139,59,80"
	read(*,*) nx,ny,nz
	dx=molxlen/(nx-1)
	dy=molylen/(ny-1)
	dz=molzlen/(nz-1)
else if (igridsel==7) then
	write(*,*) "Input X,Y,Z coordinate of original point (Bohr) e.g. 0.1,4,-1"
	read(*,*) orgx,orgy,orgz
	write(*,*) "Input X,Y,Z component of translation vector (Bohr) e.g. 0.1,0.1,0.15"
	read(*,*) dx,dy,dz
	write(*,*) "Input the number of points in X,Y,Z direction e.g. 139,59,80"
	read(*,*) nx,ny,nz
	endx=orgx+dx*(nx-1)
	endy=orgy+dy*(ny-1)
	endz=orgz+dz*(nz-1)
else if (igridsel==8) then
	write(*,*) "Input X,Y,Z coordinate of box center (in Angstrom)"
	write(*,*) "or input such as a8 to take the coordinate of atom 8 as box center"
	write(*,*) "or input such as a3,a7 to take the midpoint of atom 3 and atom 7 as box center"
	read(*,"(a)") c80tmp
	if (c80tmp(1:1)=='a') then
		do ich=1,len_trim(c80tmp)
			if (c80tmp(ich:ich)==',') exit
		end do
		if (ich==len_trim(c80tmp)+1) then
			read(c80tmp(2:),*) itmp
			cenx=a(itmp)%x
			ceny=a(itmp)%y
			cenz=a(itmp)%z
		else
			read(c80tmp(2:ich-1),*) itmp
			read(c80tmp(ich+2:),*) jtmp			
			cenx=(a(itmp)%x+a(jtmp)%x)/2D0
			ceny=(a(itmp)%y+a(jtmp)%y)/2D0
			cenz=(a(itmp)%z+a(jtmp)%z)/2D0
		end if
	else
		read(c80tmp,*) cenx,ceny,cenz
		cenx=cenx/b2a
		ceny=ceny/b2a
		cenz=cenz/b2a
	end if
	write(*,*) "Input the grid spacing (bohr)  e.g. 0.08"
	read(*,*) dx
	dy=dx
	dz=dx
	write(*,*) "Input the box lengths in X,Y,Z direction (Bohr) e.g. 8.0,8.0,13.5"
	read(*,*) molxlen,molylen,molzlen
	orgx=cenx-molxlen/2D0
	orgy=ceny-molylen/2D0
	orgz=cenz-molzlen/2D0
	endx=orgx+molxlen
	endy=orgy+molylen
	endz=orgz+molzlen
	nx=nint(molxlen/dx)+1
	ny=nint(molylen/dy)+1
	nz=nint(molzlen/dz)+1
else if (igridsel==9) then
	write(*,*) "Input path of a cube file, C:\Fate_Saber.cub"
	do while(.true.)
		read(*,"(a)") cubefilename
		inquire(file=cubefilename,exist=alive)
		if (alive) then
			open(10,file=cubefilename,status="old")
			read(10,*)
			read(10,*)
			read(10,*) nouse,orgx,orgy,orgz
			read(10,*) nx,dx
			read(10,*) ny,rnouse,dy
			read(10,*) nz,rnouse,rnouse,dz
			close(10)
			exit
		else
			write(*,*) "Error: File cannot be found, input again"
		end if
	end do
	endx=orgx+dx*(nx-1)
	endy=orgy+dy*(ny-1)
	endz=orgz+dz*(nz-1)
end if

if (igridsel==10) call setboxGUI

!If system is symmetric to Cartesian plane, then slightly adjust grid setting, so that the distribution of grids are symmetric to Cartesian plane
!This treatment will make the integrals have much better symmetry
diffx=abs(orgx+endx)
diffy=abs(orgy+endy)
diffz=abs(orgz+endz)
if (igridsel>=1.and.igridsel<=6) then
	if (diffx<0.05D0) then !The system is symmetry to YZ plane
		distxmin=1D10
		do ix=1,nx
			rnowx=orgx+(ix-1)*dx
			if (rnowx>=0D0.and.rnowx<distxmin) distxmin=rnowx
		end do
		!1D-12 is a perturbation to avoid the grids are degenerate with respect to Cartesian plane, which results in cumbersome degenerate attractors
		!However if it is introduced, the integrals will not so close to symmetry
		orgx=orgx+(dx/2D0-distxmin) !+1D-15
	end if
	if (diffy<0.05D0) then !The system is symmetry to XZ plane
		distymin=1D10
		do iy=1,ny
			rnowy=orgy+(iy-1)*dy
			if (rnowy>=0D0.and.rnowy<distymin) distymin=rnowy
		end do
		orgy=orgy+(dy/2D0-distymin) !+1D-15
	end if
	if (diffz<0.05D0) then !The system is symmetry to XY plane
		distzmin=1D10
		do iz=1,nz
			rnowz=orgz+(iz-1)*dz
			if (rnowz>=0D0.and.rnowz<distzmin) distzmin=rnowz
		end do
		orgz=orgz+(dz/2D0-distzmin) !+1D-15
	end if
end if

write(*,"(' Coordinate of origin in X,Y,Z is   ',3f12.6)") orgx,orgy,orgz
write(*,"(' Coordinate of end point in X,Y,Z is',3f12.6)") endx,endy,endz
write(*,"(' Spacing in X,Y,Z is',3f11.6)") dx,dy,dz
write(*,"(' Number of points in X,Y,Z is',3i5,'   Total',i10)") nx,ny,nz,nx*ny*nz

gridvec1=0;gridvec1(1)=dx
gridvec2=0;gridvec2(2)=dy
gridvec3=0;gridvec3(3)=dz
! do i=1,nx
! 	write(c80tmp,"(D20.13)") orgx+(i-1)*dx
! 	read(c80tmp,*) tmpval
! 	write(*,*) i," x ",tmpval
! end do
! do i=1,ny
! 	write(c80tmp,"(D20.13)") orgy+(i-1)*dy
! 	read(c80tmp,*) tmpval
! 	write(*,*) i," y ",tmpval
! end do
! do i=1,nz
! 	write(c80tmp,"(D20.13)") orgz+dfloat(i-1)*dz
! 	read(c80tmp,*) tmpval
! 	write(*,*) i," z ",tmpval
! end do
! read(*,*)
end subroutine


!!-- Detect proper box (namely set orgx/y/z,endx/y/z) so that the box can just enclose the isosurface of rho=rhocrit
subroutine detectboxrho(rhocrit)
use defvar
use function
implicit real*8(a-h,o-z)
real*8 rhocrit
tmpfac=1.0D0
orgxtmp=minval( a(:)%x-tmpfac*vdwr_tianlu(a(:)%index) ) !Define a too small box, extend each side (to obtain orgx/y/z,endx/y/z) by detecting rho
orgytmp=minval( a(:)%y-tmpfac*vdwr_tianlu(a(:)%index) )
orgztmp=minval( a(:)%z-tmpfac*vdwr_tianlu(a(:)%index) )
endxtmp=maxval( a(:)%x+tmpfac*vdwr_tianlu(a(:)%index) )
endytmp=maxval( a(:)%y+tmpfac*vdwr_tianlu(a(:)%index) )
endztmp=maxval( a(:)%z+tmpfac*vdwr_tianlu(a(:)%index) )
scanstp=0.3D0 !Spacing of scan
nstpx=nint((endxtmp-orgxtmp)/scanstp)+1 !The number of detection grid in X direction
nstpy=nint((endytmp-orgytmp)/scanstp)+1 !The number of detection grid in Y direction
nstpz=nint((endztmp-orgztmp)/scanstp)+1 !The number of detection grid in Z direction
distmove=0.2D0
write(*,*) "Detecting proper box size..."
!Detect Z low. Scan XY plane, if rho at a point is larger than "rhocrit", lower down the plane by "distmove" and check again, until all points have rho<rhocrit
orgz=orgztmp
do while(.true.)
	iok=1
zl:do ix=1,nstpx
		rnowx=orgxtmp+(ix-1)*scanstp
		do iy=1,nstpy
			rnowy=orgytmp+(iy-1)*scanstp
			if (fdens(rnowx,rnowy,orgz)>rhocrit) then
				iok=0
				exit zl
			end if
		end do
	end do zl
	if (iok==1) exit
	orgz=orgz-0.2D0 !Lower down the layer
end do
!Detect Z high
endz=endztmp
do while(.true.)
	iok=1
zh:do ix=1,nstpx
		rnowx=orgxtmp+(ix-1)*scanstp
		do iy=1,nstpy
			rnowy=orgytmp+(iy-1)*scanstp
			if (fdens(rnowx,rnowy,endz)>rhocrit) then
				iok=0
				exit zh
			end if
		end do
	end do zh
	if (iok==1) exit
	endz=endz+0.2D0
end do
!Detect X low. Scan YZ plane
orgx=orgxtmp
do while(.true.)
	iok=1
xl:do iy=1,nstpy
		rnowy=orgytmp+(iy-1)*scanstp
		do iz=1,nstpz
			rnowz=orgztmp+(iz-1)*scanstp
			if (fdens(orgx,rnowy,rnowz)>rhocrit) then
				iok=0
				exit xl
			end if
		end do
	end do xl
	if (iok==1) exit
	orgx=orgx-0.2D0
end do
!Detect X high
endx=endxtmp
do while(.true.)
	iok=1
xh:do iy=1,nstpy
		rnowy=orgytmp+(iy-1)*scanstp
		do iz=1,nstpz
			rnowz=orgztmp+(iz-1)*scanstp
			if (fdens(endx,rnowy,rnowz)>rhocrit) then
				iok=0
				exit xh
			end if
		end do
	end do xh
	if (iok==1) exit
	endx=endx+0.2D0
end do
!Detect Y low. Scan XZ plane
orgy=orgytmp
do while(.true.)
	iok=1
yl:do ix=1,nstpx
		rnowx=orgxtmp+(ix-1)*scanstp
		do iz=1,nstpz
			rnowz=orgztmp+(iz-1)*scanstp
			if (fdens(rnowx,orgy,rnowz)>rhocrit) then
				iok=0
				exit yl
			end if
		end do
	end do yl
	if (iok==1) exit
	orgy=orgy-0.2D0
end do
!Detect Y high
endy=endytmp
do while(.true.)
	iok=1
yh:do ix=1,nstpx
		rnowx=orgxtmp+(ix-1)*scanstp
		do iz=1,nstpz
			rnowz=orgztmp+(iz-1)*scanstp
			if (fdens(rnowx,endy,rnowz)>rhocrit) then
				iok=0
				exit yh
			end if
		end do
	end do yh
	if (iok==1) exit
	endy=endy+0.2D0
end do

end subroutine




!!------- Integrate a real space function in the basins already partitioned
subroutine integratebasin
use defvar
use util
use function
use basinintmod
implicit real*8(a-h,o-z)
real*8 intval(-1:numatt),basinvol(-1:numatt),intvalpriv(-1:numatt),basinvolpriv(-1:numatt)
integer walltime1,walltime2
character grdfilename*200

write(*,*) "Please select integrand:"
write(*,*) "-2 Return"
write(*,*) "-1 The values of the grid data stored in an external file (.cub/.grd)"
write(*,*) "0 The values of the grid data stored in memory"
call selfunc_interface(1,ifuncint)
if (ifuncint==-1) then
	do while(.true.)
		write(*,*) "Input another .cub or .grd file name, e.g. C:\CODE_GEASS\CC.cub"
		read(*,"(a)") grdfilename
		inquire(file=grdfilename,exist=alive)
		if (alive) exit
		write(*,*) "File not found, input again"
		write(*,*)
	end do
	inamelen=len_trim(grdfilename)
	if (grdfilename(inamelen-2:inamelen)=="cub".or.grdfilename(inamelen-3:inamelen)=="cube") then
		call readcubetmp(grdfilename,1,inconsis)
	else if (grdfilename(inamelen-2:inamelen)=="grd") then
		call readgrdtmp(grdfilename,inconsis)
	end if
	if (inconsis==1) return
	write(*,*)
else if (ifuncint==-2) then
	return
end if

call walltime(walltime1)
write(*,*) "Integrating, please wait..."
intval=0D0
basinvol=0D0
ifinish=0
!$OMP PARALLEL private(ix,iy,iz,irealatt,rnowx,rnowy,rnowz,tmpval,intvalpriv,basinvolpriv) shared(intval,basinvol,ifinish) NUM_THREADS(nthreads)
intvalpriv=0D0
basinvolpriv=0D0
!$OMP do schedule(DYNAMIC)
do iz=2,nz-1
	do iy=2,ny-1
		do ix=2,nx-1
			rnowx=orgx+(ix-1)*dx
			rnowy=orgy+(iy-1)*dy
			rnowz=orgz+(iz-1)*dz
			if (ifuncint==-1) then
				tmpval=cubmattmp(ix,iy,iz)
			else if (ifuncint==0) then
				tmpval=cubmat(ix,iy,iz)
			else
				tmpval=calcfuncall(ifuncint,rnowx,rnowy,rnowz)
			end if
			irealatt=gridbas(ix,iy,iz)
			intvalpriv(irealatt)=intvalpriv(irealatt)+tmpval
			basinvolpriv(irealatt)=basinvolpriv(irealatt)+1
		end do
	end do
    ifinish=ifinish+1
    if (ifuncint>=1) call showprog(ifinish,nz-2)
end do
!$OMP end do
!$OMP CRITICAL
    intval=intval+intvalpriv
    basinvol=basinvol+basinvolpriv
!$OMP end CRITICAL
!$OMP END PARALLEL
dvol=dx*dy*dz
intval=intval*dvol
basinvol=basinvol*dvol !Basin volume
write(*,*) "  #Basin        Integral(a.u.)      Volume(a.u.^3)"
do irealatt=1,numrealatt
	write(*,"(i8,f22.10,f20.8)") irealatt,intval(irealatt),basinvol(irealatt)
end do
write(*,"(' Sum of above values:',f20.8)") sum(intval(1:numrealatt))
if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==0)) write(*,"(' Integral of unassigned grids:',f20.8)") intval(0)
if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==-1)) write(*,"(' Integral of the grids travelled to box boundary:',f20.8)") intval(-1)

call walltime(walltime2)
write(*,"(' Integrating basins took up wall clock time',i10,' s')") walltime2-walltime1
end subroutine



!!------- Calculate orbital compositions contributed by various basins by uniform grid integration
subroutine basinorbcomp
use defvar
use util
use basinintmod
use function
implicit real*8(a-h,o-z)
real*8 intval(-1:numatt),intvalpriv(-1:numatt),orbcomp(ncenter)
real*8 comparr(nmo,numrealatt),comparrpriv(nmo,numrealatt),orbval(nmo)
integer att2atm(numrealatt),atm2att(ncenter)
character c80tmp*80,c2000tmp*2000

if (allocated(frag1)) deallocate(frag1)
if (ifuncbasin==1) then
    call calc_att2atm(att2atm)
    atm2att=0 !If an atom doesn't have corresponding attractor, the attractor index will be 0
    do irealatt=1,numrealatt
        atm2att(att2atm(irealatt))=irealatt
    end do
end if
dvol=dx*dy*dz

do while(.true.)
    write(*,*)
	write(*,*) "Now input the orbital index to print orbital composition, e.g. 5"
	write(*,"(a)") " You can also input:"
	if (.not.allocated(frag1)) then
	    write(*,"(a,i6)") " -9: Define fragment, current: undefined"
	else
    	write(*,"(a,i6)") " -9: Redefine fragment, current number of atoms:",size(frag1)
    end if
	write(*,"(a)") "  0: Return"
	write(*,"(a)") " -1: Print basic information of all orbitals"
	!write(*,"(a)") " -2: Print atom contribution to a batch of orbitals"
	!write(*,"(a)") " -3: Print fragment contribution to a batch of orbitals"
	write(*,"(a)") " -4: Export composition of every basin in every orbital to orbcomp.txt"
	!write(*,"(a)") " -5: Print orbital delocalization index (ODI) for a batch of orbitals"
    
	read(*,*) isel
    if (isel==0) then
        if (allocated(frag1)) deallocate(frag1)
        return
    else if (isel==-1) then
        call showorbinfo(1,nmo)
        
    else if (isel==-4) then
        if (ifuncbasin/=1) then
            write(*,"(a)") " Error: This option is meaningful only when the function used to partition the basins is electron density!"
            write(*,*) "Press ENTER button to continue"
            read(*,*)
            cycle
        end if
        comparr=0D0
        ifinish=0
        !$OMP PARALLEL private(ix,iy,iz,irealatt,rnowx,rnowy,rnowz,comparrpriv,orbval) shared(ifinish,comparr) NUM_THREADS(nthreads)
        comparrpriv=0D0
        !$OMP do schedule(DYNAMIC)
        do iz=2,nz-1
	        do iy=2,ny-1
		        do ix=2,nx-1
			        rnowx=orgx+(ix-1)*dx
			        rnowy=orgy+(iy-1)*dy
			        rnowz=orgz+(iz-1)*dz
				    call orbderv(1,1,nmo,rnowx,rnowy,rnowz,orbval(:))
			        irealatt=gridbas(ix,iy,iz)
			        comparrpriv(:,irealatt)=comparrpriv(:,irealatt)+orbval(:)**2
		        end do
	        end do
            ifinish=ifinish+1
            call showprog(ifinish,nz-2)
        end do
        !$OMP end do
        !$OMP CRITICAL
            comparr=comparr+comparrpriv
        !$OMP end CRITICAL
        !$OMP END PARALLEL
        comparr=comparr*dvol
    
		open(10,file="orbcomp.txt",status="replace")
		do imo=1,nmo
			write(10,"(' Orbital',i6)") imo
            do iatm=1,ncenter
                orbcomp(iatm)=comparr(imo,atm2att(iatm))
            end do
            valnorm=sum(comparr(imo,:)) !Since there may be NNA, normalization factor should use sum of all basins
			do iatm=1,ncenter
				write(10,"(i6,f11.3,' %')") iatm,orbcomp(iatm)/valnorm*100
			end do
		end do
		close(10)
		write(*,*) "Done! orbcomp.txt has been exported to current folder"
        if (numrealatt/=ncenter) then
            write(*,"(a)") " Note: The number of attractors is unequal to number of atoms, therefore the sum of atom contributions is not equal to 100 %"
        end if
        
    else if (isel==-9) then
        if (ifuncbasin/=1) then
            write(*,"(a)") " Error: This option is meaningful only when the function used to partition the basins is electron density!"
            write(*,*) "Press ENTER button to continue"
            read(*,*)
            cycle
        end if
        call definefragment
        
    else
        intval=0D0
        ifinish=0
        !$OMP PARALLEL private(ix,iy,iz,irealatt,rnowx,rnowy,rnowz,tmpval,intvalpriv) shared(intval,ifinish) NUM_THREADS(nthreads)
        intvalpriv=0D0
        !$OMP do schedule(DYNAMIC)
        do iz=2,nz-1
	        do iy=2,ny-1
		        do ix=2,nx-1
			        rnowx=orgx+(ix-1)*dx
			        rnowy=orgy+(iy-1)*dy
			        rnowz=orgz+(iz-1)*dz
				    tmpval=fmo(rnowx,rnowy,rnowz,isel)**2
			        irealatt=gridbas(ix,iy,iz)
			        intvalpriv(irealatt)=intvalpriv(irealatt)+tmpval
		        end do
	        end do
            ifinish=ifinish+1
            call showprog(ifinish,nz-2)
        end do
        !$OMP end do
        !$OMP CRITICAL
            intval=intval+intvalpriv
        !$OMP end CRITICAL
        !$OMP END PARALLEL
        intval=intval*dvol
        sumcontri=sum(intval(1:numrealatt))
        write(*,"(' Sum of raw contributions:',f18.3,' %')") sumcontri*100
        intval=intval/sumcontri
        write(*,*)
        write(*,*) "Final data after normalization:"
        do irealatt=1,numrealatt
	        write(*,"(' Basin:',i8,'   Contribution:',f10.3,' %')") irealatt,intval(irealatt)*100
        end do
        !For AIM basins, also show atomic values
        if (ifuncbasin==1) then
            write(*,*)
            write(*,*) "Contributions from atoms:"
            do iatm=0,ncenter
			    do iatt=1,numrealatt
				    if (att2atm(iatt)==iatm) then
                        if (iatm==0) then
					        write(*,"(i7,' (NNA)   Contribution:',f8.3,' %')") iatt,intval(iatt)*100
				        else
				            write(*,"(i7,' (',a,')    Contribution:',f8.3,' %')") iatm,a(iatm)%name,intval(iatt)*100
				        end if
                    end if
			    end do
		    end do
            write(*,*)
            write(*,"(' Orbital delocalization index:',f8.2)") sum((intval(1:numrealatt)*100)**2)/100
		    if (allocated(frag1)) then
			    write(*,*)
                compfrag=sum(intval(atm2att(frag1)))
			    write(*,"(' Composition of the fragment:',f8.3,' %')") compfrag*100
                orbdeloc=sum((intval(atm2att(frag1))*100/compfrag)**2)
                write(*,"(' Orbital delocalization index of the fragment:',f8.2)") orbdeloc/100
		    end if
        end if
    end if
end do
end subroutine




!!------- Obtain atomic population in a basin region (e.g. ELF bond basin) via AIM partition
subroutine atmpopinbasin
use defvar
use basinintmod
implicit real*8(a-h,o-z)
inquire(file="basin.cub",exist=alive)
if (.not.alive) then
	write(*,*) "Error: basin.cub is not existed in current folder!"
	return
else
	if (allocated(cubmattmp)) deallocate(cubmattmp)
	call readcubetmp("basin.cub",1,inconsis)
	if (inconsis==1) then
		write(*,*) "Error: The grid setting of basin.cub is inconsistent with present grid data"
		return
	end if
end if
do while(.true.)
	write(*,*)
	write(*,*) "Study population of which atom? Input the index of corresponding attractor"
	write(*,*) "For example, 3"
	write(*,*) "Input 0 can exit"
	read(*,*) iatt
	if (iatt==0) then
		exit
	else if (iatt<0.or.iatm>numrealatt) then
		write(*,*) "Error: Attractor index exceeded valid range!"
		return
	end if
	write(*,*) "Study its population in which basin of the basin.cub file? e.g. 6"
	read(*,*) ibasin
	atmpop=0
	do iz=2,nz-1
		do iy=2,ny-1
			do ix=2,nx-1
				if (nint(cubmattmp(ix,iy,iz))==ibasin.and.gridbas(ix,iy,iz)==iatt) atmpop=atmpop+cubmat(ix,iy,iz)
			end do
		end do
	end do
	atmpop=atmpop*dx*dy*dz
	write(*,"(' Population of attractor',i4,' in external basin',i4,' is',f12.5)") iatt,ibasin,atmpop
end do
deallocate(cubmattmp)
end subroutine



!!----------------- Calculate multipole moment in the basins
subroutine multipolebasin
use defvar
use basinintmod
use function
implicit real*8 (a-h,o-z)
do while(.true.)
	write(*,*) "Input the index of the basin, e.g. 5"
	write(*,"(a)") " Note: -1 means printing all basin results on screen, -2 means printing to multipol.txt in current folder. Input 0 can return"
	read(*,*) itmp
	if (itmp==-1.or.itmp==-2) then
		ibegin=1
		iend=numrealatt
		if (itmp==-2) then
			ioutid=10
			open(10,file="multipol.txt",status="replace")
		end if
	else if (itmp==0) then
		return
	else
		ibegin=itmp
		iend=itmp
		ioutid=6
	end if
	write(*,*) "Calculating, please wait..."
	write(ioutid,*) "Note: All units shown below are in a.u.!"
	write(ioutid,*)
	
	dipelextot=0D0
	dipeleytot=0D0
	dipeleztot=0D0
	eleinttot=0D0
	do ibas=ibegin,iend
	! 	xcen=0 !Use centroid of basin as center, but this is a bad idea
	! 	ycen=0
	! 	zcen=0
	! 	nbasgrid=0
	! 	do iz=2,nz-1
	! 		do iy=2,ny-1
	! 			do ix=2,nx-1
	! 				if (gridbas(ix,iy,iz)==ibas) then
	! 					nbasgrid=nbasgrid+1
	! 					rnowx=orgx+(ix-1)*dx
	! 					rnowy=orgy+(iy-1)*dy
	! 					rnowz=orgz+(iz-1)*dz
	! 					tmpdens=ELF_LOL(rnowx,rnowy,rnowz,"ELF")
	! 					xcen=xcen+rnowx*tmpdens
	! 					ycen=ycen+rnowy*tmpdens
	! 					zcen=zcen+rnowz*tmpdens
	! 				end if
	! 			end do
	! 		end do
	! 	end do
	! 	xcen=xcen/nbasgrid
	! 	ycen=ycen/nbasgrid
	! 	zcen=zcen/nbasgrid
		xcen=realattxyz(ibas,1)
		ycen=realattxyz(ibas,2)
		zcen=realattxyz(ibas,3)
		dvol=dx*dy*dz
	! 	write(*,"(' The X,Y,Z of the center of the basin:')")
	! 	write(*,"(3f14.8,' Bohr',/)") xcen,ycen,zcen
	
		eleint=0D0
		xint=0D0
		yint=0D0
		zint=0D0
		xxint=0D0
		yyint=0D0
		zzint=0D0
		xyint=0D0
		yzint=0D0
		xzint=0D0
		!$OMP PARALLEL private(ix,iy,iz,rnowx,rnowy,rnowz,tmpmul,rx,ry,rz,eleintp,xintp,yintp,zintp,xxintp,yyintp,zzintp,xyintp,yzintp,xzintp) NUM_THREADS(nthreads)
		eleintp=0D0
		xintp=0D0
		yintp=0D0
		zintp=0D0
		xxintp=0D0
		yyintp=0D0
		zzintp=0D0
		xyintp=0D0
		yzintp=0D0
		xzintp=0D0
		!$OMP do schedule(DYNAMIC)
		do iz=2,nz-1
			rnowz=orgz+(iz-1)*dz
			do iy=2,ny-1
				rnowy=orgy+(iy-1)*dy
				do ix=2,nx-1
					if (gridbas(ix,iy,iz)==ibas) then
						rnowx=orgx+(ix-1)*dx
						if (ifuncbasin==1) then
							tmpmul=dvol*cubmat(ix,iy,iz) !The cubmat currently is just electron density
						else
							tmpmul=dvol*fdens(rnowx,rnowy,rnowz)
						end if
						rx=rnowx-xcen
						ry=rnowy-ycen
						rz=rnowz-zcen
						eleintp=eleintp+tmpmul !monopole
						xintp=xintp+rx*tmpmul
						yintp=yintp+ry*tmpmul
						zintp=zintp+rz*tmpmul
						xxintp=xxintp+rx*rx*tmpmul
						yyintp=yyintp+ry*ry*tmpmul
						zzintp=zzintp+rz*rz*tmpmul
						xyintp=xyintp+rx*ry*tmpmul
						yzintp=yzintp+ry*rz*tmpmul
						xzintp=xzintp+rx*rz*tmpmul
					end if
				end do
			end do
		end do
		!$OMP end do
		!$OMP CRITICAL
		eleint=eleint+eleintp
		xint=xint+xintp
		yint=yint+yintp
		zint=zint+zintp
		xxint=xxint+xxintp
		yyint=yyint+yyintp
		zzint=zzint+zzintp
		xyint=xyint+xyintp
		yzint=yzint+yzintp
		xzint=xzint+xzintP
		!$OMP end CRITICAL
		!$OMP END PARALLEL
		if (itmp==-1.or.itmp==-2) write(ioutid,"(' Basin',i8)") ibas
		if (itmp==-2) write(*,"(' Outputting basin',i8)") ibas
		write(ioutid,"(' Basin electric monopole moment:',f12.6)") -eleint
		write(ioutid,"(' Basin electric dipole moment:')") 
		write(ioutid,"('  X=',f12.6,'  Y=',f12.6,'  Z=',f12.6,'  Magnitude=',f12.6)") -xint,-yint,-zint,sqrt(xint**2+yint**2+zint**2)
		eleinttot=eleinttot+eleint
		dipelex=-eleint*xcen+(-xint) !Contribution to molecular total dipole moment
		dipeley=-eleint*ycen+(-yint)
		dipelez=-eleint*zcen+(-zint)
		dipelextot=dipelextot+dipelex
		dipeleytot=dipeleytot+dipeley
		dipeleztot=dipeleztot+dipelez
		write(ioutid,"(' Basin electron contribution to molecular dipole moment:')")
		write(ioutid,"('  X=',f12.6,'  Y=',f12.6,'  Z=',f12.6,'  Magnitude=',f12.6)") dipelex,dipeley,dipelez,sqrt(dipelex**2+dipeley**2+dipelez**2)
		write(ioutid,"(' Basin electric quadrupole moment (Cartesian form):')")
		rrint=xxint+yyint+zzint
		QXX=-(3*xxint-rrint)/2
		QYY=-(3*yyint-rrint)/2
		QZZ=-(3*zzint-rrint)/2
		write(ioutid,"(' QXX=',f12.6,'  QXY=',f12.6,'  QXZ=',f12.6)") QXX,-(3*xyint)/2,-(3*xzint)/2
		write(ioutid,"(' QYX=',f12.6,'  QYY=',f12.6,'  QYZ=',f12.6)") -(3*xyint)/2,QYY,-(3*yzint)/2
		write(ioutid,"(' QZX=',f12.6,'  QZY=',f12.6,'  QZZ=',f12.6)") -(3*xzint)/2,-(3*yzint)/2,QZZ
		write(ioutid,"( ' The magnitude of electric quadrupole moment (Cartesian form):',f12.6)") dsqrt(2D0/3D0*(QXX**2+QYY**2+QZZ**2))
		R20=-(3*zzint-rrint)/2D0 !Notice that the negative sign, because electrons carry negative charge
		R2n1=-dsqrt(3D0)*yzint
		R2p1=-dsqrt(3D0)*xzint
		R2n2=-dsqrt(3D0)*xyint
		R2p2=-dsqrt(3D0)/2D0*(xxint-yyint)
		write(ioutid,"(' Electric quadrupole moments (Spherical harmonic form):')")
		write(ioutid,"(' Q_2,0 =',f11.6,'   Q_2,-1=',f11.6,'   Q_2,1=',f11.6)") R20,R2n1,R2p1
		write(ioutid,"(' Q_2,-2=',f11.6,'   Q_2,2 =',f11.6)") R2n2,R2p2
		write(ioutid,"( ' Magnitude: |Q_2|=',f12.6)") dsqrt(R20**2+R2n1**2+R2p1**2+R2n2**2+R2p2**2)
		write(ioutid,*)
	end do
	if (itmp==-1.or.itmp==-2) then !Output overall properties, most users are not interested in them
		dipnucx=sum(a(:)%x*a(:)%charge)
		dipnucy=sum(a(:)%y*a(:)%charge)
		dipnucz=sum(a(:)%z*a(:)%charge)
		write(ioutid,"( ' Molecular net charge:',f12.6)") sum(a%charge)-eleinttot
		write(ioutid,"( ' Nuclear contribution to molecular dipole moment:')") 
		write(ioutid,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") dipnucx,dipnucy,dipnucz,sqrt(dipnucx**2+dipnucy**2+dipnucz**2)
		write(ioutid,"( ' Electron contribution to molecular dipole moment:')") 
		write(ioutid,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") dipelextot,dipeleytot,dipeleztot,sqrt(dipelextot**2+dipeleytot**2+dipeleztot**2)
		dipmolx=dipnucx+dipelextot
		dipmoly=dipnucy+dipeleytot
		dipmolz=dipnucz+dipeleztot
		write(ioutid,"( ' Molecular dipole moment:')")
		write(ioutid,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") dipmolx,dipmoly,dipmolz,sqrt(dipmolx**2+dipmoly**2+dipmolz**2)
		write(ioutid,*)
	end if
	if (itmp==-2) then
		close(10)
		write(*,*) "Outputting finished!"
		write(*,*)
	end if
end do
end subroutine


!!------------------ Calculate localized index and delocalization index within and between basins
!Most of the codes are adapted from fuzzyana.f90
!itask==0 means calculate BOM and LI,DI
!itask==1 means calculate and output BOM to BOM.txt in current folder
!itask==2 means calculate and output AOM to AOM.txt in current folder
!BOM means overlap matrix of all occupied orbitals in basins
subroutine LIDIbasin(itask)
use defvar
use basinintmod
use function
use util
implicit real*8 (a-h,o-z)
real*8 DI(numrealatt,numrealatt),DIa(numrealatt,numrealatt),DIb(numrealatt,numrealatt) !Delocalization index matrix
real*8 LI(numrealatt),LIa(numrealatt),LIb(numrealatt) !Localization index array
real*8,allocatable :: BOM(:,:,:),BOMa(:,:,:),BOMb(:,:,:),BOMsum(:,:),BOMsuma(:,:),BOMsumb(:,:) !BOM(i,j,k) means overlap matrix of MO i,j in atom k space
real*8 :: BOMtmp(nmo,nmo),orbval(nmo),atmDIa(ncenter,ncenter),atmDIb(ncenter,ncenter)
character selectyn
integer itask,maplist(ncenter)

!Allocate space for BOM
if (wfntype==0.or.wfntype==2.or.wfntype==3) then !RHF,ROHF,R-post-HF
	if (wfntype==3) then !R-post-HF, need to consider all orbitals
		nmatsize=nmo
	else !RHF,ROHF
		!High-lying virtual orbitals will be deleted, especially for .fch case (In fact, before entering this function, those >= LUMO+10 has already been discarded)
		!Notice that occupation number may be not contiguous, some low-lying orbitals may have
		!zero occupation due to modification by users, so we can't simply use nelec to determine matrix size
		do nmatsize=nmo,1,-1
			if (MOocc(nmatsize)/=0) exit
		end do
		if (nmo-nmatsize>0) write(*,"(' Note: The highest',i6,' virtual orbitals will not be taken into account')") nmo-nmatsize
	end if
	if (allocated(BOM)) deallocate(BOM,BOMsum)
	allocate(BOM(nmatsize,nmatsize,numrealatt),BOMsum(nmatsize,nmatsize))
	BOM=0
	BOMsum=0
else if (wfntype==1.or.wfntype==4) then !UHF,U-post-HF
	do iendalpha=nmo,1,-1
		if (MOtype(iendalpha)==1) exit
	end do
	if (wfntype==4) then
		nmatsizea=iendalpha !Total number of alpha orbitals
		nmatsizeb=nmo-nmatsizea !Total number of beta orbitals
	else
		do nmatsizea=iendalpha,1,-1
			if (MOocc(nmatsizea)/=0D0) exit
		end do
		if (nint(nbelec)==0) then
			nmatsizeb=0
		else
			do nmatsizeb=nmo,iendalpha+1,-1
				if (MOocc(nmatsizeb)/=0D0) exit
			end do
			nmatsizeb=nmatsizeb-iendalpha
		end if
		if (iendalpha-nmatsizea>0) write(*,"(' Note: The highest',i6,' alpha virtual orbitals will not be taken into account')") iendalpha-nmatsizea
		if (nmo-iendalpha-nmatsizeb>0) write(*,"(' Note: The highest',i6,' beta virtual orbitals will not be taken into account')") nmo-iendalpha-nmatsizeb
	end if
	if (allocated(BOMa)) deallocate(BOMa,BOMb,BOMsuma,BOMsumb)
	allocate( BOMa(nmatsizea,nmatsizea,numrealatt),BOMb(nmatsizeb,nmatsizeb,numrealatt) )
	allocate( BOMsuma(nmatsizea,nmatsizea),BOMsumb(nmatsizeb,nmatsizeb) )
	BOMa=0
	BOMb=0
	BOMsuma=0
	BOMsumb=0
end if

!Calculate BOM
dvol=dx*dy*dz
call walltime(nwalltime1)
if (wfntype==0.or.wfntype==2.or.wfntype==3) then !RHF,ROHF,R-post-HF
	do ibas=1,numrealatt
		write(*,"(' Generating orbital overlap matrix for basin',i6,'  of',i6,' ......')") ibas,numrealatt
		!$OMP parallel shared(BOM) private(ix,iy,iz,rnowx,rnowy,rnowz,imo,jmo,BOMtmp,orbval) num_threads(nthreads)
		BOMtmp=0D0
		!$OMP do schedule(DYNAMIC)
		do iz=2,nz-1
			do iy=2,ny-1
				do ix=2,nx-1
					if (gridbas(ix,iy,iz)==ibas) then
						rnowx=orgx+(ix-1)*dx
						rnowy=orgy+(iy-1)*dy
						rnowz=orgz+(iz-1)*dz
						call orbderv(1,1,nmatsize,rnowx,rnowy,rnowz,orbval)
						do imo=1,nmatsize
							do jmo=imo,nmatsize
								BOMtmp(imo,jmo)=BOMtmp(imo,jmo)+orbval(imo)*orbval(jmo)*dvol
							end do
						end do
					end if
				end do
			end do
		end do
		!$OMP end do
		!$OMP CRITICAL
		BOM(:,:,ibas)=BOM(:,:,ibas)+BOMtmp(1:nmatsize,1:nmatsize)
		!$OMP end CRITICAL
		!$OMP end parallel
		BOM(:,:,ibas)=BOM(:,:,ibas)+transpose(BOM(:,:,ibas))
		do imo=1,nmatsize
			BOM(imo,imo,ibas)=BOM(imo,imo,ibas)/2D0
		end do
	end do
else if (wfntype==1.or.wfntype==4) then !UHF,U-post-HF
	do ibas=1,numrealatt
		!Alpha part
		write(*,"(' Generating orbital overlap matrix for basin',i6,'  of',i6,' ......')") ibas,numrealatt
		!$OMP parallel shared(BOMa) private(ix,iy,iz,rnowx,rnowy,rnowz,imo,jmo,BOMtmp,orbval) num_threads(nthreads)
		BOMtmp=0D0
		!$OMP do schedule(DYNAMIC)
		do iz=2,nz-1
			do iy=2,ny-1
				do ix=2,nx-1
					if (gridbas(ix,iy,iz)==ibas) then
						rnowx=orgx+(ix-1)*dx
						rnowy=orgy+(iy-1)*dy
						rnowz=orgz+(iz-1)*dz
						call orbderv(1,1,nmatsizea,rnowx,rnowy,rnowz,orbval)
						do imo=1,nmatsizea
							do jmo=imo,nmatsizea
								BOMtmp(imo,jmo)=BOMtmp(imo,jmo)+orbval(imo)*orbval(jmo)*dvol
							end do
						end do
					end if
				end do
			end do
		end do
		!$OMP end do
		!$OMP CRITICAL
		BOMa(:,:,ibas)=BOMa(:,:,ibas)+BOMtmp(1:nmatsizea,1:nmatsizea)
		!$OMP end CRITICAL
		!$OMP end parallel
		BOMa(:,:,ibas)=BOMa(:,:,ibas)+transpose(BOMa(:,:,ibas))
		do imo=1,nmatsizea
			BOMa(imo,imo,ibas)=BOMa(imo,imo,ibas)/2D0
		end do
		
		!Beta part
		if (nmatsizeb>0) then !If there is no beta orbital, then nmatsizeb=0, and we will do nothing
			MOinit=iendalpha+1
			MOend=iendalpha+nmatsizeb
! 			write(*,*) MOinit,MOend,nmatsizeb
			!$OMP parallel shared(BOMb) private(ix,iy,iz,rnowx,rnowy,rnowz,imo,imotmp,jmo,jmotmp,BOMtmp,orbval) num_threads(nthreads)
			BOMtmp=0D0
			!$OMP do schedule(DYNAMIC)
			do iz=2,nz-1
				do iy=2,ny-1
					do ix=2,nx-1
						if (gridbas(ix,iy,iz)==ibas) then
							rnowx=orgx+(ix-1)*dx
							rnowy=orgy+(iy-1)*dy
							rnowz=orgz+(iz-1)*dz
							call orbderv(1,MOinit,MOend,rnowx,rnowy,rnowz,orbval)
							do imo=MOinit,MOend
								imotmp=imo-iendalpha !So that the index start from 1 to nbelec
								do jmo=imo,MOend
									jmotmp=jmo-iendalpha
									BOMtmp(imotmp,jmotmp)=BOMtmp(imotmp,jmotmp)+orbval(imo)*orbval(jmo)*dvol
								end do
							end do
						end if
					end do
				end do
			end do
			!$OMP end do
			!$OMP CRITICAL
			BOMb(:,:,ibas)=BOMb(:,:,ibas)+BOMtmp(1:nmatsizeb,1:nmatsizeb)
			!$OMP end CRITICAL
			!$OMP end parallel
			BOMb(:,:,ibas)=BOMb(:,:,ibas)+transpose(BOMb(:,:,ibas))
			do imo=1,nmatsizeb
				BOMb(imo,imo,ibas)=BOMb(imo,imo,ibas)/2D0
			end do
		end if
	end do
end if

write(*,*)
!Check sanity of BOM
if (wfntype==0.or.wfntype==2.or.wfntype==3) then !RHF,ROHF,R-post-HF
	do ibas=1,numrealatt
		BOMsum=BOMsum+BOM(:,:,ibas)
	end do
	BOMerror=identmaterr(BOMsum)/numrealatt
	write(*,"(' Error of BOM is',f14.8)") BOMerror
	if (BOMerror>0.02D0) write(*,"(a)") " Warning: The integration is not very accurate"
! 	call showmatgau(BOMsum,"BOMsum",0,"f14.8",6)
else if (wfntype==1.or.wfntype==4) then !UHF,U-post-HF
	BOMerrorb=0D0
	do ibas=1,numrealatt
		BOMsuma=BOMsuma+BOMa(:,:,ibas)
		if (nmatsizeb>0) BOMsumb=BOMsumb+BOMb(:,:,ibas)
	end do
	BOMerrora=identmaterr(BOMsuma)/numrealatt
	if (nmatsizeb>0) BOMerrorb=identmaterr(BOMsumb)/numrealatt
	write(*,"(' Error of alpha BOM is',f14.8)") BOMerrora
	if (nmatsizeb>0) write(*,"(' Error of Beta BOM is ',f14.8)") BOMerrorb
	if (BOMerrora>0.02D0.or.BOMerrorb>0.02D0) write(*,"(a)") " Warning: The integration is not very accurate"
! 	call showmatgau(BOMsumb,"BOMsumb",0,"f14.8",6)
end if

if (itask==1) then !Output BOM
	open(10,file="BOM.txt",status="replace")
	if (wfntype==0.or.wfntype==2.or.wfntype==3) then
		do ibas=1,numrealatt
			write(10,"('Orbital overlap matrix of basin',i6)") ibas
			call showmatgau(BOM(:,:,ibas),"",1,"f14.8",10)
			write(10,*)
		end do
	else if (wfntype==1.or.wfntype==4) then
		do ibas=1,numrealatt
			write(10,"('Alpha part of orbital overlap matrix of basin',i6)") ibas
			call showmatgau(BOMa(:,:,ibas),"",1,"f14.8",10)
			if (nmatsizeb>0) then
				write(10,"('Beta part of orbital overlap matrix of basin',i6)") ibas
				call showmatgau(BOMb(:,:,ibas),"",1,"f14.8",10)
			end if
			write(10,*)
		end do
	end if
	close(10)
	write(*,*)
	write(*,*) "Done, the matrices have been exported to BOM.txt in current folder"
	return
else if (itask==2) then !Output AOM
    write(*,*)
    call atmidx2attidx(maplist,1,6)
	open(10,file="AOM.txt",status="replace")
    do iatm=1,ncenter
        ibas=maplist(iatm)
	    if (wfntype==0.or.wfntype==2.or.wfntype==3) then
		    write(10,"('Atomic overlap matrix of',i6,'(',a,')')") iatm,a(iatm)%name
		    call showmatgau(BOM(:,:,ibas),"",1,"f14.8",10)
		    write(10,*)
	    else if (wfntype==1.or.wfntype==4) then
		    write(10,"('Alpha part of atomic overlap matrix of',i6,'(',a,')')") iatm,a(iatm)%name
		    call showmatgau(BOMa(:,:,ibas),"",1,"f14.8",10)
		    if (nmatsizeb>0) then
			    write(10,"('Beta part of atomic overlap matrix of',i6,'(',a,')')") iatm,a(iatm)%name
			    call showmatgau(BOMb(:,:,ibas),"",1,"f14.8",10)
		    end if
		    write(10,*)
        end if
    end do
	close(10)
	write(*,*)
	write(*,*) "Done, the matrices have been exported to AOM.txt in current folder"
	return
end if

!Generate DI and LI
!RHF,R-post-HF, DI_A,B=2¡Æ[i,j]dsqrt(n_i*n_j)*S_i,j_A * S_i,j_B     where i and j are non-spin orbitals
write(*,*) "Generating LI and DI..."
if (any(MOocc<0)) then
	where(MOocc<0) MOocc=0
	write(*,"(a)") " Note: Some occupation numbers are negative. In order to make the calculation feasible, they have been set to zero"
	write(*,*) "Press ENTER button to continue"
	read(*,*)
end if
if (wfntype==0.or.wfntype==3) then
	DI=0D0
	do ibas=1,numrealatt
		do jbas=ibas,numrealatt
			do iorb=1,nmatsize
				do jorb=1,nmatsize
					DI(ibas,jbas)=DI(ibas,jbas)+dsqrt(MOocc(iorb)*MOocc(jorb))*BOM(iorb,jorb,ibas)*BOM(iorb,jorb,jbas)
				end do
			end do
		end do
		LI(ibas)=DI(ibas,ibas)
	end do
	DI=2*(DI+transpose(DI))
	do ibas=1,numrealatt !Diagonal terms are the sum of corresponding row or column
		DI(ibas,ibas)=0D0
		DI(ibas,ibas)=sum(DI(ibas,:))
	end do
else if (wfntype==2) then !ROHF
	DIa=0D0
	DIb=0D0
	do nmoclose=nmatsize,1,-1
		if (MOtype(nmoclose)==0) exit
	end do
	do ibas=1,numrealatt
		do jbas=ibas,numrealatt
			!Alpha
			do iorb=1,nmatsize !The number of close or alpha orbitals needed to be concerned
				occi=MOocc(iorb)
				if (MOtype(iorb)==0) occi=occi/2D0
				do jorb=1,nmatsize
					occj=MOocc(jorb)
					if (MOtype(jorb)==0) occj=occj/2D0
					DIa(ibas,jbas)=DIa(ibas,jbas)+dsqrt(occi*occj)*BOM(iorb,jorb,ibas)*BOM(iorb,jorb,jbas)
				end do
			end do
			!Beta
			do iorb=1,nmoclose !The number of close orbitals needed to be concerned
				do jorb=1,nmoclose
					DIb(ibas,jbas)=DIb(ibas,jbas)+dsqrt(MOocc(iorb)/2D0*MOocc(jorb)/2D0)*BOM(iorb,jorb,ibas)*BOM(iorb,jorb,jbas)
				end do
			end do
		end do
		LIa(ibas)=DIa(ibas,ibas)
		LIb(ibas)=DIb(ibas,ibas)
	end do
	DIa=2*(DIa+transpose(DIa))
	DIb=2*(DIb+transpose(DIb))
	do ibas=1,numrealatt !Diagonal terms are the sum of corresponding row or column
		DIa(ibas,ibas)=0D0
		DIb(ibas,ibas)=0D0
		DIa(ibas,ibas)=sum(DIa(ibas,:))
		DIb(ibas,ibas)=sum(DIb(ibas,:))
	end do
	!Combine alpha and Beta to total
	DI=DIa+DIb
	LI=LIa+LIb
!UHF,U-post-HF   DI(A,B)=2¡Æ[i,j]dsqrt(n_i*n_j)*S_i,j_A * S_i,j_B   where i and j are spin orbitals
else if (wfntype==1.or.wfntype==4) then
	!Alpha
	DIa=0D0
	do ibas=1,numrealatt
		do jbas=ibas,numrealatt
			do iorb=1,nmatsizea
				do jorb=1,nmatsizea
					DIa(ibas,jbas)=DIa(ibas,jbas)+dsqrt(MOocc(iorb)*MOocc(jorb))*BOMa(iorb,jorb,ibas)*BOMa(iorb,jorb,jbas)
				end do
			end do
		end do
		LIa(ibas)=DIa(ibas,ibas)
	end do
	DIa=2*(DIa+transpose(DIa))
	!Beta
	if (nmatsizeb>0) then
		DIb=0D0
		MOinit=iendalpha+1 !Index range of beta orbitals
		MOend=iendalpha+nmatsizeb
		do ibas=1,numrealatt
			do jbas=ibas,numrealatt
				do iorb=MOinit,MOend
					iorbtmp=iorb-iendalpha
					do jorb=MOinit,MOend
						jorbtmp=jorb-iendalpha
						DIb(ibas,jbas)=DIb(ibas,jbas)+dsqrt(MOocc(iorb)*MOocc(jorb))*BOMb(iorbtmp,jorbtmp,ibas)*BOMb(iorbtmp,jorbtmp,jbas)
					end do
				end do
			end do
			LIb(ibas)=DIb(ibas,ibas)
		end do
		DIb=2*(DIb+transpose(DIb))
	end if
	do ibas=1,numrealatt !Diagonal terms are the sum of corresponding row or column
		DIa(ibas,ibas)=0D0
		DIb(ibas,ibas)=0D0
		DIa(ibas,ibas)=sum(DIa(ibas,:))
		DIb(ibas,ibas)=sum(DIb(ibas,:))
	end do
	!Combine alpha and Beta to total
	DI=DIa+DIb
	LI=LIa+LIb
end if

call walltime(nwalltime2)
write(*,"(' Calculation took up',i8,' seconds wall clock time',/)")  nwalltime2-nwalltime1

!Output LI and DI
ioutid=6
100 write(ioutid,"(a,/)") " Note: Diagonal terms are the sum of corresponding row or column elements"
if (ifuncbasin==1) write(ioutid,"(a,/)") " Note: In below output the indices are basin indices rather than atom indices"
if (wfntype==1.or.wfntype==2.or.wfntype==4) then !UHF,ROHF,U-post-HF, output each spin component first
	!Alpha
	call showmatgau(DIa,"Delocalization index matrix for alpha spin",0,"f14.8",ioutid)
	write(ioutid,*)
	write(ioutid,*) "Localization index for alpha electron:"
	do ibas=1,numrealatt
		write(ioutid,"(i6,':',f7.3)",advance='no') ibas,LIa(ibas)
		if (mod(ibas,5)==0) write(ioutid,*)
	end do
    if (mod(numrealatt,5)/=0) write(ioutid,*)
	write(ioutid,*)
	!Beta
	call showmatgau(DIb,"Delocalization index matrix for beta spin",0,"f14.8",ioutid)
	write(ioutid,*)
	write(ioutid,*) "Localization index for beta spin:"
	do ibas=1,numrealatt
		write(ioutid,"(i6,':',f7.3)",advance='no') ibas,LIb(ibas)
		if (mod(ibas,5)==0) write(ioutid,*)
	end do
    if (mod(numrealatt,5)/=0) write(ioutid,*)
	write(ioutid,*)
end if
!Alpha+Beta
call showmatgau(DI,"Total delocalization index matrix",0,"f14.8",ioutid)
write(ioutid,*)
write(ioutid,*) "Total localization index:"
do ibas=1,numrealatt
	write(ioutid,"(i6,':',f7.3)",advance='no') ibas,LI(ibas)
	if (mod(ibas,5)==0) write(ioutid,*)
end do
if (mod(numrealatt,5)/=0) write(ioutid,*)

!Write LI and DI in atom indices
if (ifuncbasin==1) then
    write(ioutid,*)
    call atmidx2attidx(maplist,1,ioutid)
    write(ioutid,"(/,a,/)") " Note: In below output the indices are atom indices rather than basin indices"
    if (wfntype==1.or.wfntype==2.or.wfntype==4) then !UHF,ROHF,U-post-HF, output each spin component first
	    !Alpha
        atmDIa=0
        atmDIb=0
        do iatm=1,ncenter
            do jatm=1,ncenter
                if (maplist(iatm)==0.or.maplist(jatm)==0) cycle !Corresponding basin index was not identified
                atmDIa(iatm,jatm)=DIa(maplist(iatm),maplist(jatm))
                atmDIb(iatm,jatm)=DIb(maplist(iatm),maplist(jatm))
            end do
        end do
	    call showmatgau(atmDIa,"Delocalization index matrix (atom index) for alpha spin",0,"f14.4",ioutid,formindex="6x,i5,3x")
	    write(ioutid,*)
	    write(ioutid,*) "Localization index (atom index) for alpha electron:"
	    do iatm=1,ncenter
            if (maplist(iatm)==0) cycle
		    write(ioutid,"(i6,':',f7.3)",advance='no') iatm,LIa(maplist(iatm))
		    if (mod(iatm,5)==0) write(ioutid,*)
	    end do
	    write(ioutid,*)
	    write(ioutid,*)
	    !Beta
	    call showmatgau(atmDIb,"Delocalization index matrix (atom index) for beta spin",0,"f14.4",ioutid,formindex="6x,i5,3x")
	    write(ioutid,*)
	    write(ioutid,*) "Localization index (atom index) for beta spin:"
	    do iatm=1,ncenter
            if (maplist(iatm)==0) cycle
		    write(ioutid,"(i6,':',f7.3)",advance='no') ibas,LIb(maplist(iatm))
		    if (mod(iatm,5)==0) write(ioutid,*)
	    end do
	    write(ioutid,*)
	    write(ioutid,*)
    end if
    !Alpha+Beta
    atmDIa=0
    do iatm=1,ncenter
        do jatm=1,ncenter
            if (maplist(iatm)==0.or.maplist(jatm)==0) cycle !Corresponding basin index was not identified
            atmDIa(iatm,jatm)=DI(maplist(iatm),maplist(jatm))
        end do
    end do
    call showmatgau(atmDIa,"Total delocalization index matrix (atom index)",0,"f14.4",ioutid,formindex="6x,i5,3x")
    write(ioutid,*)
    write(ioutid,*) "Total localization index (atom index):"
    do iatm=1,ncenter
        if (maplist(iatm)==0) cycle
	    write(ioutid,"(i6,':',f7.3)",advance='no') iatm,LI(maplist(iatm))
	    if (mod(iatm,5)==0) write(ioutid,*)
    end do
    if (mod(numrealatt,5)/=0) write(ioutid,*)
end if

if (ioutid==10) then
	write(*,*) "Done!"
	close(10)
	return
end if

write(*,*)
write(*,*) "If also outputting above information to LIDI.txt in current folder? (y/n)"
read(*,*) selectyn
if (selectyn=='y'.or.selectyn=='Y') then
	open(10,file="LIDI.txt",status="replace")
	ioutid=10
	goto 100
end if
end subroutine



!!------- Find basin indices corresponding to various atom indices in the case of AIM analysis
!maplist(i) is the attractor index that corresponds to atom i. If no corresponding can be found, the value is 0
!info: =1 show correspondence, =0 do not show
!ioutid: Destination of output information
!Distance <0.3 Bohr is the criterion used to identify correspondence
subroutine atmidx2attidx(maplist,info,ioutid)
use defvar
use basinintmod
implicit real*8 (a-h,o-z)
integer maplist(ncenter),ioutid
distcrit=0.3D0
maplist=0
if (info==1) write(ioutid,*) "Detecting correspondence between basin and atom indices (criterion: <0.3 Bohr)"
do iatm=1,ncenter
	do ibas=1,numrealatt
        disttest=dsqrt( (realattxyz(ibas,1)-a(iatm)%x)**2+(realattxyz(ibas,2)-a(iatm)%y)**2+(realattxyz(ibas,3)-a(iatm)%z)**2 )
        if (disttest<distcrit) then
            if (info==1) write(ioutid,"(' Basin',i6,' corresponds to atom',i6,' (',a,')')") ibas,iatm,a(iatm)%name 
            maplist(iatm)=ibas
            exit
        end if
    end do
    if (ibas==numrealatt+1.and.info==1) write(ioutid,"(a,i6)") " Warning: Unable to find basin index that corresponds to atom",iatm
end do
end subroutine



!!------- Integrate a real space function in the basins with multi-level refinement, indenspensible for Laplacian
!DEPRECATED, since mixed type of grids perform better and faster
subroutine integratebasinrefine
use defvar
use util
use function
use basinintmod
implicit real*8(a-h,o-z)
character c200tmp*200
real*8 intval(-1:numatt),basinvol(-1:numatt),intvalpriv(-1:numatt),basinvolpriv(-1:numatt)
integer walltime1,walltime2
real*8 :: critlevel1=0.1D0,critlevel2=0.5D0,critlevel3=1D0
integer :: nrefine1=1,nrefine2=2,nrefine3=5,nrefine4=7
if (ifuncbasin/=1) then
	write(*,*) "Error: This function is only applicable to AIM basins!"
	return
end if

do while(.true.)
	write(*,*) "-2 Return"
	write(*,*) "-1 Print and set parameters for multi-level refinement"
	call selfunc_interface(1,ifuncint)
	if (ifuncint==-2) then
		return
	else if (ifuncint==-1) then
		write(*,"(a)") " Note: A number n means each grid in corresponding value range will be transformed &
		as n^3 grids around it during the integration to gain a higher integration accuracy. n=1 means the grids will remain unchanged. &
		The ""value range"" referred here is the value range of the function used to generate basins"
		write(*,*) "The value range and the times of refinement:"
		write(*,"(' Smaller than ',f10.5,' :',i4)") critlevel1,nrefine1
		write(*,"(' Between',f10.5,' and ',f10.5,' :',i4)") critlevel1,critlevel2,nrefine2
		write(*,"(' Between',f10.5,' and ',f10.5,' :',i4)") critlevel2,critlevel3,nrefine3
		write(*,"(' Larger than  ',f10.5,' :',i4)") critlevel3,nrefine4
		write(*,*)
		write(*,*) "Please input three thresholds to define the four ranges, e.g. 0.1,0.5,1"
		write(*,*) "Note: Press ENTER button can retain current values unchanged"
		read(*,"(a)") c200tmp
		if (c200tmp/=' ') read(c200tmp,*) critlevel1,critlevel2,critlevel3
		write(*,*) "Input the times of refinement for the grids in the four ranges, e.g. 1,2,5,7"
		write(*,*) "Note: Press ENTER button can retain current values unchanged"
		read(*,"(a)") c200tmp
		if (c200tmp/=' ') read(c200tmp,*) nrefine1,nrefine2,nrefine3,nrefine4
		write(*,*) "Done!"
		write(*,*)
	end if
end do

call walltime(walltime1)
write(*,*) "Integrating, please wait..."
intval=0D0
basinvol=0D0
ifinish=0
!$OMP PARALLEL private(ix,iy,iz,ixref,iyref,izref,ndiv,irealatt,rnowx,rnowy,rnowz,rnowxtmp,rnowytmp,rnowztmp,orgxref,orgyref,orgzref,dxref,dyref,dzref,&
!$OMP tmpval,tmpvalrefine,intvalpriv,basinvolpriv,nrefine) shared(intval,basinvol,ifinish) NUM_THREADS(nthreads)
intvalpriv=0D0
basinvolpriv=0D0
!$OMP do schedule(DYNAMIC)
do iz=2,nz-1
	do iy=2,ny-1
		do ix=2,nx-1
			if (cubmat(ix,iy,iz)<critlevel1) then
				nrefine=nrefine1 !The number of point to represent each edge
			else if (cubmat(ix,iy,iz)<critlevel2) then
				nrefine=nrefine2
			else if (cubmat(ix,iy,iz)<critlevel3) then
				nrefine=nrefine3
			else
				nrefine=nrefine4
			end if
			ndiv=nrefine**3
			rnowx=orgx+(ix-1)*dx
			rnowy=orgy+(iy-1)*dy
			rnowz=orgz+(iz-1)*dz
			orgxref=rnowx-dx/2 !Take corner position as original point of microcycle
			orgyref=rnowy-dy/2
			orgzref=rnowz-dz/2
			dxref=dx/nrefine
			dyref=dy/nrefine
			dzref=dz/nrefine
			tmpval=0D0
			do ixref=1,nrefine
				do iyref=1,nrefine
					do izref=1,nrefine
						rnowxtmp=orgxref+(ixref-0.5D0)*dxref
						rnowytmp=orgyref+(iyref-0.5D0)*dyref
						rnowztmp=orgzref+(izref-0.5D0)*dzref
						if (ifuncint==-1) then
							tmpvalrefine=cubmattmp(ix,iy,iz)
						else if (ifuncint==0) then
							tmpvalrefine=cubmat(ix,iy,iz)
						else
							tmpvalrefine=calcfuncall(ifuncint,rnowxtmp,rnowytmp,rnowztmp)
						end if
						tmpval=tmpval+tmpvalrefine/ndiv
					end do
				end do
			end do
			irealatt=gridbas(ix,iy,iz)
			intvalpriv(irealatt)=intvalpriv(irealatt)+tmpval
			basinvolpriv(irealatt)=basinvolpriv(irealatt)+1
		end do
	end do
    ifinish=ifinish+1
    call showprog(ifinish,nz-2)
end do
!$OMP end do
!$OMP CRITICAL
    intval=intval+intvalpriv
    basinvol=basinvol+basinvolpriv
!$OMP end CRITICAL
!$OMP END PARALLEL
dvol=dx*dy*dz
intval=intval*dvol
basinvol=basinvol*dvol !Basin volume
write(*,*) "  #Basin          Integral        Volume(a.u.^3)"
do irealatt=1,numrealatt
	write(*,"(i8,f22.10,f20.8)") irealatt,intval(irealatt),basinvol(irealatt)
end do
write(*,"(' Sum of above values:',f20.8)") sum(intval(1:numrealatt))
if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==0)) write(*,"(' Integral of unassigned grids:',f20.8)") intval(0)
if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==-1)) write(*,"(' Integral of the grids travelled to box boundary:',f20.8)") intval(-1)

call walltime(walltime2)
write(*,"(' Integrating basins took up wall clock time',i10,' s')") walltime2-walltime1
end subroutine



!!------- Integrate AIM basins using mixed atomic-center and uniform grids
! itype=1: Integrate specific real space function
! itype=2: Integrate specific real space function with exact refinement of basin boundary
! itype=3: Integrate specific real space function with approximate refinement of basin boundary by exact+linear interpolation
! itype=10: Produce electric multipole moments
!NNA and ECP are supported
!Notice that the grid generated must cover the whole space!
subroutine integratebasinmix(itype)
use defvar
use util
use function
use basinintmod
use topo
implicit real*8(a-h,o-z)
!p suffix means private variable for parallel mode
real*8 intval(-1:numrealatt,20),intvalp(-1:numrealatt,20) !Up to 20 functions can be evaluated and stored simultaneously
real*8 basinvol(-1:numrealatt),basinvolp(-1:numrealatt),basinvdwvol(-1:numrealatt),basinvdwvolp(-1:numrealatt) !vdW is used to obtain the basin volume enclosed by 0.001 isosurface of rho
real*8 trustrad(numrealatt),intbasinthread(numrealatt),intbasin(numrealatt)
real*8 dens,grad(3),hess(3,3),k1(3),k2(3),k3(3),k4(3),xarr(nx),yarr(ny),zarr(nz)
real*8,allocatable :: potx(:),poty(:),potz(:),potw(:)
real*8,allocatable :: rhogrid(:,:,:),rhogradgrid(:,:,:,:) !Used in shubin's 2nd project
real*8,allocatable :: prorhogrid(:,:,:) !Used in integrating deformation density
type(content),allocatable :: gridatt(:) !Record correspondence between attractor and grid
integer att2atm(numrealatt) !The attractor corresponds to which atom. If =0, means this is a NNA
real*8 eleint(-1:numrealatt),xint(-1:numrealatt),yint(-1:numrealatt),zint(-1:numrealatt),&
xxint(-1:numrealatt),yyint(-1:numrealatt),zzint(-1:numrealatt),xyint(-1:numrealatt),yzint(-1:numrealatt),xzint(-1:numrealatt)
real*8 eleintp(-1:numrealatt),xintp(-1:numrealatt),yintp(-1:numrealatt),zintp(-1:numrealatt),&
xxintp(-1:numrealatt),yyintp(-1:numrealatt),zzintp(-1:numrealatt),xyintp(-1:numrealatt),yzintp(-1:numrealatt),xzintp(-1:numrealatt)
integer walltime1,walltime2,radpotAIM,sphpotAIM
real*8 prodensgrad(0:4),gridval(100000,6) !Used for storing various information during integrating inside sphere
character c80tmp*80,selectyn
nbeckeiter=8
if (ifuncbasin/=1) then
	write(*,"(a)") " Error: This function is only applicable to AIM basins! That means in option 1, you should select electron density to construct the basins."
	return
end if

if (ispecial==0) then
	if (itype==1.or.itype==2.or.itype==3) then
		write(*,*) "Please select integrand:"
		write(*,*) "-2 Return"
		write(*,*) "-1 Deformation density"
		call selfunc_interface(1,ifuncint)
		if (ifuncint==-2) then
			return
		else if (ifuncint==-1) then
			call setpromol
		end if
	end if
else if (ispecial==1) then
	continue !Don't let user to select integrand
else if (ispecial==2) then
	call setpromol !Don't let user to select integrand
	expcutoff=1 !Use full accuracy for shubin's 2nd project
end if

call walltime(walltime1)
CALL CPU_TIME(time_begin)
intval=0D0
basinvol=0D0
basinvdwvol=0D0
eleint=0D0
xint=0D0
yint=0D0
zint=0D0
xxint=0D0
yyint=0D0
zzint=0D0
xyint=0D0
yzint=0D0
xzint=0D0
numcp=0

att2atm=0
!Determine trust radius and then integrate in the trust sphere
!We use CPpos to record the center of the trust sphere
write(*,*) "Integrating in trust sphere..."
do iatt=1,numrealatt !Cycle each attractors
	!Refine the raw position of attractor by exact Newton method
	do iatm=1,ncenter
		disttest=dsqrt( (realattxyz(iatt,1)-a(iatm)%x)**2+(realattxyz(iatt,2)-a(iatm)%y)**2+(realattxyz(iatt,3)-a(iatm)%z)**2 )
		if (disttest<0.3D0) then !If distance between attractor and a nucleus is smaller than 0.3 Bohr, then the attractor will belong to the atom
			att2atm(iatt)=iatm
			write(*,"(' Attractor',i6,' corresponds to atom',i6,' (',a,')')") iatt,iatm,a(iatm)%name
			numcpold=numcp
			call findcp(a(iatm)%x,a(iatm)%y,a(iatm)%z,1)
			if (numcp==numcpold) then
				write(*,*) "Note: Unable to locate exact CP position! Use nuclear position"
				numcp=numcp+1
				CPpos(1,numcp)=a(iatm)%x
				CPpos(2,numcp)=a(iatm)%y
				CPpos(3,numcp)=a(iatm)%z
			end if
! 			write(*,"(' Coordinate after refinement:',3f16.8)") CPpos(:,numcp)
			exit
		end if
	end do
	if (att2atm(iatt)==0) then
		write(*,"(a,i6,a)") " Warning: Unable to determine the attractor",iatt," belongs to which atom!"
		write(*,"(a)") " If this is a non-nuclear attractor, simply press ENTER button to continue. If you used pseudopotential &
		and this attractor corresponds to the cluster of all maxima of its valence electron, then input the index of this atom (e.g. 9). &
		Else you should input q to return and regenerate basins with smaller grid spacing"
		read(*,"(a)") c80tmp
		if (c80tmp=='q') then
			return
		else if (c80tmp==" ") then
			numcpold=numcp
			call findcp(realattxyz(iatt,1),realattxyz(iatt,2),realattxyz(iatt,3),1)
			if (numcp==numcpold) then
				write(*,*) "Unable to locate exact CP position! Exit..."
				return
			end if
		else !ECP, input the corresponding atom by user
			read(c80tmp,*) iatmtmp
			att2atm(iatt)=iatmtmp
			numcp=numcp+1
			CPpos(1,numcp)=a(iatmtmp)%x
			CPpos(2,numcp)=a(iatmtmp)%y
			CPpos(3,numcp)=a(iatmtmp)%z
		end if
	end if
	
	!Determine trust radius and set integration points and weight
	radpotAIM=200
	parm=1
	isettrustrad=0
	nintgrid=0 !Then number of integration grids within trust radius
	if (allocated(gridatt)) deallocate(gridatt) !Used to record grids in trust sphere of this attractor
	allocate(gridatt(radpotAIM*500))
	do ish=1,radpotAIM !Cycle each radial shell. Radius distance is from near to far
		if (isettrustrad==1) exit !The trust radius has been finally determined in last shell cycle
		!Becke, namely the second-kind Gauss-Chebyshev
		itmp=radpotAIM+1-ish !Invert ish to make radr from near to far
		radx=cos(itmp*pi/(radpotAIM+1D0))
		radr=(1+radx)/(1-radx)*parm
		radw=2*pi/(radpotAIM+1)*parm**3 *(1+radx)**2.5D0/(1-radx)**3.5D0 *4*pi
! 		!Handy, also known as Euler-Maclaurin. See: Murray, C. W.; Handy, N. C.; Laming, G. J. Mol Phys 1993, 78, 997
! 		radx=dfloat(ish)/(radpotAIM+1D0)
! 		radr=radx**2/(1-radx)**2*parm
! 		radw=2*radx**5/dfloat(radpotAIM+1)/(1-radx)**7*parm**3 *4*pi
		
		!Set Lebedev grids according to shell radius
		!For more inner shell, the distribution is more akin to spherically symmetric, therefore lower number of grids could be used
		if (att2atm(iatt)==0) then !NNA
			sphpotAIM=302
		else
			radtmp=covr(a(att2atm(iatt))%index)
			if (radr<0.2D0*radtmp) then
				sphpotAIM=26
			else if (radr<0.5D0*radtmp) then
				sphpotAIM=74
			else if (radr<0.8D0*radtmp) then
				sphpotAIM=146
			else
				sphpotAIM=194
			end if
		end if
		if (allocated(potx)) deallocate(potx,poty,potz,potw)
		allocate(potx(sphpotAIM),poty(sphpotAIM),potz(sphpotAIM),potw(sphpotAIM))
		call Lebedevgen(sphpotAIM,potx,poty,potz,potw)
		!Combine radial point and weights with angular part, and make them centered at current attractor
		gridatt( nintgrid+1:nintgrid+sphpotAIM )%x=radr*potx+CPpos(1,numcp)
		gridatt( nintgrid+1:nintgrid+sphpotAIM )%y=radr*poty+CPpos(2,numcp)
		gridatt( nintgrid+1:nintgrid+sphpotAIM )%z=radr*potz+CPpos(3,numcp)
		gridatt( nintgrid+1:nintgrid+sphpotAIM )%value=radw*potw
		!Find out trust radius for present attractor
		!If in a shell, the angle between "linking line between nucleus and a shell point" and "gradient vector of this point" &
		!is larger than 45 degree, then this shell is trust radius
		angmax=0
		if (att2atm(iatt)==0) then
			radinit=0
		else
			radrinit=0.15D0
			if (a(att2atm(iatt))%index>2) radrinit=0.5D0
		end if
		if (isettrustrad==0.and.radr>radrinit) then
			do isphpt=1,sphpotAIM
				xtmp=gridatt(nintgrid+isphpt)%x
				ytmp=gridatt(nintgrid+isphpt)%y
				ztmp=gridatt(nintgrid+isphpt)%z
				call calchessmat_dens(1,xtmp,ytmp,ztmp,dens,grad,hess) !Only value and gradient
				dirx=CPpos(1,numcp)-xtmp
				diry=CPpos(2,numcp)-ytmp
				dirz=CPpos(3,numcp)-ztmp
				angtmp=vecang(dirx,diry,dirz,grad(1),grad(2),grad(3))
				if (angtmp>angmax) angmax=angtmp
				if (angtmp>45) then
					write(*,"(' The trust radius of attractor',i6,' is',f10.3,' Bohr',/)") iatt,trustrad(iatt)
					isettrustrad=1 !The radius of last shell should be the final trust radius. Now exit
					exit
				end if
			end do
			if (isettrustrad==0) trustrad(iatt)=radr !Passed this shell and temporarily set the radius as trust radius. Continue to enlarge the trust radius, until reached angmax>45 degree
		end if
		nintgrid=nintgrid+sphpotAIM
	end do
	if (isettrustrad==0) then !Trust radius was not set after run over all shells
		trustrad(iatt)=1000 !Infinite, for isolated atom
		if (ispecial==2) trustrad(iatt)=20 !For Shubin's 2nd project, should not be as large as 1000, because for a point very far from nucleus the relative entropy cannot be evaluated
		write(*,"(' The trust radius of attractor',i6,' is',f10.3,' Bohr',/)") iatt,trustrad(iatt)
	end if
	
	!Use DFT integration algorithm to integrate the region inside trust radius
	!$OMP PARALLEL private(ipt,ptx,pty,ptz,rx,ry,rz,dist,tmps,iter,switchwei,intvalp,&
	!$OMP eleintp,xintp,yintp,zintp,xxintp,yyintp,zzintp,xyintp,yzintp,xzintp,tmpval,tmpval2,tmpval3,tmpval4) &
	!$OMP shared(intval,gridval) NUM_THREADS(nthreads)
	intvalp=0D0
	eleintp=0D0
	xintp=0D0
	yintp=0D0
	zintp=0D0
	xxintp=0D0
	yyintp=0D0
	zzintp=0D0
	xyintp=0D0
	yzintp=0D0
	xzintp=0D0
	!$OMP do schedule(DYNAMIC)
	do ipt=1,nintgrid
		ptx=gridatt(ipt)%x
		pty=gridatt(ipt)%y
		ptz=gridatt(ipt)%z
		rx=ptx-CPpos(1,numcp) !The relative distance between current point to corresponding attractor
		ry=pty-CPpos(2,numcp)
		rz=ptz-CPpos(3,numcp)
		!Calculate switching function
		dist=dsqrt(rx*rx+ry*ry+rz*rz)
		tmps=dist-trustrad(iatt)
		if (tmps>1) then
			switchwei=0
		else if (tmps<-1) then
			switchwei=1
		else
			do iter=1,nbeckeiter
				tmps=1.5D0*(tmps)-0.5D0*(tmps)**3
			end do
			switchwei=0.5D0*(1-tmps)
		end if
		gridval(ipt,3)=switchwei
! 		if (dist>trustrad(iatt)) cycle !Discrete separation between atomic-center and uniform integration
		if (switchwei<1D-7) cycle !For saving computational time
		
		if (itype==1.or.itype==2.or.itype==3) then !Integrate a function
			if (ifuncint==-1) then !Deformation density, store molecular density of present center temporarily
				gridval(ipt,1)=fdens(ptx,pty,ptz)
			else if (ispecial==0) then !Normal case
				tmpval=calcfuncall(ifuncint,ptx,pty,ptz)
				intvalp(iatt,1)=intvalp(iatt,1)+gridatt(ipt)%value*tmpval*switchwei
			else if (ispecial==1) then !For Shubin's project, simultaneously output many properties
				tmpval=infoentro(2,ptx,pty,ptz) !Shannon entropy density
				tmpval2=Fisherinfo(1,ptx,pty,ptz) !Fisher information density
				tmpval3=weizsacker(ptx,pty,ptz) !Steric energy
				tmpval4=fdens(ptx,pty,ptz) !Electron density
				intvalp(iatt,1)=intvalp(iatt,1)+gridatt(ipt)%value*tmpval*switchwei
				intvalp(iatt,2)=intvalp(iatt,2)+gridatt(ipt)%value*tmpval2*switchwei
				intvalp(iatt,3)=intvalp(iatt,3)+gridatt(ipt)%value*tmpval3*switchwei
				intvalp(iatt,4)=intvalp(iatt,4)+gridatt(ipt)%value*tmpval4*switchwei
			else if (ispecial==2) then !For Shubin's 2nd project
				call calchessmat_dens(1,ptx,pty,ptz,gridval(ipt,1),gridval(ipt,4:6),hess)
				gridval(ipt,2)=sum(gridval(ipt,4:6)**2)
			end if
		else if (itype==10) then !Calculate multipole moment
			tmpval=gridatt(ipt)%value*fdens(ptx,pty,ptz)*switchwei
			eleintp(iatt)=eleintp(iatt)+tmpval
			xintp(iatt)=xintp(iatt)+rx*tmpval
			yintp(iatt)=yintp(iatt)+ry*tmpval
			zintp(iatt)=zintp(iatt)+rz*tmpval
			xxintp(iatt)=xxintp(iatt)+rx*rx*tmpval
			yyintp(iatt)=yyintp(iatt)+ry*ry*tmpval
			zzintp(iatt)=zzintp(iatt)+rz*rz*tmpval
			xyintp(iatt)=xyintp(iatt)+rx*ry*tmpval
			yzintp(iatt)=yzintp(iatt)+ry*rz*tmpval
			xzintp(iatt)=xzintp(iatt)+rx*rz*tmpval
		end if
	end do
	!$OMP end do
	!$OMP CRITICAL
	if (itype==1.or.itype==2.or.itype==3) then
		intval=intval+intvalp
	else if (itype==10) then
		eleint=eleint+eleintp
		xint=xint+xintp
		yint=yint+yintp
		zint=zint+zintp
		xxint=xxint+xxintp
		yyint=yyint+yyintp
		zzint=zzint+zzintp
		xyint=xyint+xyintp
		yzint=yzint+yzintp
		xzint=xzint+xzintp
	end if
	!$OMP end CRITICAL
	!$OMP END PARALLEL
	
	!Some special cases:
	if (ispecial==2) then !Shubin's 2nd project, integrate relative Shannon and Fisher information
		!Current gridval content: 1=rho, 2=gradrho^2, 3=switchwei, 4=gradrho_x, 5=gradrho_y, 6=gradrho_z
		call dealloall
		write(*,"(' Loading ',a,/)") trim(custommapname(att2atm(iatt)))
		call readwfn(custommapname(att2atm(iatt)),1)
		do ipt=1,nintgrid
			switchwei=gridval(ipt,3)
			if (switchwei<1D-7) cycle
			ptx=gridatt(ipt)%x
			pty=gridatt(ipt)%y
			ptz=gridatt(ipt)%z
			call calchessmat_dens(1,ptx,pty,ptz,prodensgrad(0),prodensgrad(1:3),hess)
			prodens=prodensgrad(0) !rho0_A
			prodensgrad2=sum(prodensgrad(1:3)**2)
			!Relative Shannon entropy, integrate rho*log(rho/rho0_A)]
			tmpval=gridval(ipt,1)*log(gridval(ipt,1)/prodens)
			intval(iatt,1)=intval(iatt,1)+gridatt(ipt)%value*tmpval*switchwei
			!Relative Fisher information (old and incorrect implementation)
			tmpval=gridval(ipt,2)/gridval(ipt,1)-prodensgrad2/prodens
			intval(iatt,2)=intval(iatt,2)+gridatt(ipt)%value*tmpval*switchwei
			!Relative Fisher information (new and correct implementation)
			tmpvalx=gridval(ipt,4)/gridval(ipt,1) - prodensgrad(1)/prodens
			tmpvaly=gridval(ipt,5)/gridval(ipt,1) - prodensgrad(2)/prodens
			tmpvalz=gridval(ipt,6)/gridval(ipt,1) - prodensgrad(3)/prodens
			tmpval=gridval(ipt,1)*(tmpvalx**2+tmpvaly**2+tmpvalz**2)
			intval(iatt,3)=intval(iatt,3)+gridatt(ipt)%value*tmpval*switchwei
		end do
		call dealloall
		call readinfile(firstfilename,1) !Retrieve to first loaded file(whole molecule) to calc real rho again
	else if (ifuncint==-1) then !Integrate deformation density
		gridval(:,2)=0D0
		do iatm=1,ncenter_org !Cycle each atom to calculate deformation density at all integration grid
			call dealloall
			call readwfn(custommapname(iatm),1)
			do ipt=1,nintgrid
				switchwei=gridval(ipt,3)
				if (switchwei<1D-7) cycle
				gridval(ipt,2)=gridval(ipt,2)+fdens(gridatt(ipt)%x,gridatt(ipt)%y,gridatt(ipt)%z)
			end do
		end do
		do ipt=1,nintgrid !Now gridval(ipt,1/2/3) records actual, promolecular density and weight at ipt
			defdens=gridval(ipt,1)-gridval(ipt,2)
			intval(iatt,1)=intval(iatt,1)+gridatt(ipt)%value*defdens*gridval(ipt,3)
		end do
		call dealloall
		call readinfile(firstfilename,1) !Retrieve to first loaded file(whole molecule) to calc real rho again
	end if
end do !End cycle attractors

if (itype==1.or.itype==2.or.itype==3) then
	write(*,*) "Integration result inside trust spheres"
	if (ispecial/=2) then
		write(*,*) "  #Sphere       Integral(a.u.)"
		do iatt=1,numrealatt
			write(*,"(i8,f22.10)") iatt,intval(iatt,1)
		end do
		write(*,"(' Sum of above values:',f20.8)") sum(intval(1:numrealatt,1))
	else if (ispecial==2) then !Shubin's 2nd project
		write(*,*) "  #Sphere      Rel.Shannon         Rel.Fisher(old)     Rel.Fisher(new)"
		do iatt=1,numrealatt
			write(*,"(i8,3f20.10)") iatt,intval(iatt,1:3)
		end do
		write(*,"(' Sum of rel.Shannon:',f20.8)") sum(intval(1:numrealatt,1))
		write(*,"(' Sum of rel.Fisher(old): ',f20.8)") sum(intval(1:numrealatt,2))
		write(*,"(' Sum of rel.Fisher(new): ',f20.8)") sum(intval(1:numrealatt,3))
	end if
end if

!Set coordinate of uniform grids
dvol=dx*dy*dz
do ix=1,nx
	xarr(ix)=orgx+(ix-1)*dx
end do
do iy=1,ny
	yarr(iy)=orgy+(iy-1)*dy
end do
do iz=1,nz
	zarr(iz)=orgz+(iz-1)*dz
end do

!--------- Integrating uniform grids, basin boundary grids will be calculated in the later stage
!NOTE: For shubin's 2nd project or deformation density, do not integrate here. At next stage boundary grids will be updated, and then at the next stage,&
!they will be integrated by a special module. Because at each grid if we reload atomic .wfn file will be too time consuming
if (ispecial==2.or.ifuncint==-1) goto 10 
write(*,*)
write(*,*) "Integrating uniform grids..."
ifinish=0
!$OMP PARALLEL private(ix,iy,iz,iatt,icp,rnowx,rnowy,rnowz,rx,ry,rz,tmpval,tmpval2,tmpval3,intvalp,basinvolp,basinvdwvolp,dist,tmps,iter,switchwei,&
!$OMP eleintp,xintp,yintp,zintp,xxintp,yyintp,zzintp,xyintp,xzintp,yzintp) shared(intval,basinvol,basinvdwvol,ifinish) NUM_THREADS(nthreads)
intvalp=0D0
basinvolp=0D0
basinvdwvolp=0D0
eleintp=0D0
xintp=0D0
yintp=0D0
zintp=0D0
xxintp=0D0
yyintp=0D0
zzintp=0D0
xyintp=0D0
yzintp=0D0
xzintp=0D0
!$OMP do schedule(DYNAMIC)
do iz=2,nz-1
	rnowz=zarr(iz)
	do iy=2,ny-1
		rnowy=yarr(iy)
		do ix=2,nx-1
			rnowx=xarr(ix)
			if ((itype==2.or.itype==3).and.interbasgrid(ix,iy,iz)) cycle !If refine boundary grid at next stage, we don't calculate them at present stage
! 			if (iz==nint(nz/2D0)) write(1,"('C',3f14.8)") rnowx*b2a,rnowy*b2a,rnowz*b2a !Examine grid distribution
			iatt=gridbas(ix,iy,iz)
! 			do icp=1,numcp
! 				dist=dsqrt( (rnowx-CPpos(1,icp))**2+(rnowy-CPpos(2,icp))**2+(rnowz-CPpos(3,icp))**2 )
! 				if (disttest<trustrad(icp)) cycle cycix !The function inside trust radius is integrated by DFT integration
! 			end do
			rx=rnowx-CPpos(1,iatt) !The relative distance between current point to corresponding attractor
			ry=rnowy-CPpos(2,iatt)
			rz=rnowz-CPpos(3,iatt)
			!Calculate switching function at current grid
			dist=dsqrt(rx*rx+ry*ry+rz*rz)
			tmps=dist-trustrad(iatt)
			if (tmps>1) then
				switchwei=0
			else if (tmps<-1) then
				switchwei=1
			else
				do iter=1,nbeckeiter
					tmps=1.5D0*(tmps)-0.5D0*(tmps)**3
				end do
				switchwei=0.5D0*(1-tmps)
			end if
			switchwei=1-switchwei
			basinvolp(iatt)=basinvolp(iatt)+1 !Calculate basin volume
			if (cubmat(ix,iy,iz)>0.001D0) basinvdwvolp(iatt)=basinvdwvolp(iatt)+1
			if (switchwei<1D-7) cycle !For saving time
			if (itype==1.or.itype==2.or.itype==3) then
				if (ispecial==0) then
					if (ifuncint==1) then !Electron density on each grid has already been calculated
						tmpval=cubmat(ix,iy,iz)
					else
						tmpval=calcfuncall(ifuncint,rnowx,rnowy,rnowz)
					end if
					intvalp(iatt,1)=intvalp(iatt,1)+tmpval*switchwei
				else if (ispecial==1) then !For RCY
					tmpval=infoentro(2,rnowx,rnowy,rnowz) !Shannon entropy density, see JCP,126,191107 for example
					tmpval2=Fisherinfo(1,rnowx,rnowy,rnowz) !Fisher information density, see JCP,126,191107 for example
					tmpval3=weizsacker(rnowx,rnowy,rnowz) !Steric energy
					tmpval4=fdens(rnowx,rnowy,rnowz) !Electron density
					intvalp(iatt,1)=intvalp(iatt,1)+tmpval*switchwei
					intvalp(iatt,2)=intvalp(iatt,2)+tmpval2*switchwei
					intvalp(iatt,3)=intvalp(iatt,3)+tmpval3*switchwei
					intvalp(iatt,4)=intvalp(iatt,4)+tmpval4*switchwei
				end if
			else if (itype==10) then
				tmpval=cubmat(ix,iy,iz)*switchwei
				eleintp(iatt)=eleintp(iatt)+tmpval
				xintp(iatt)=xintp(iatt)+rx*tmpval
				yintp(iatt)=yintp(iatt)+ry*tmpval
				zintp(iatt)=zintp(iatt)+rz*tmpval
				xxintp(iatt)=xxintp(iatt)+rx*rx*tmpval
				yyintp(iatt)=yyintp(iatt)+ry*ry*tmpval
				zzintp(iatt)=zzintp(iatt)+rz*rz*tmpval
				xyintp(iatt)=xyintp(iatt)+rx*ry*tmpval
				yzintp(iatt)=yzintp(iatt)+ry*rz*tmpval
				xzintp(iatt)=xzintp(iatt)+rx*rz*tmpval
			end if
		end do
	end do
    ifinish=ifinish+1
    if (ifuncint/=1) call showprog(ifinish,nz-2)
end do
!$OMP end do
!$OMP CRITICAL
if (itype==1.or.itype==2.or.itype==3) then
	intval=intval+intvalp*dvol
	basinvol=basinvol+basinvolp*dvol
	basinvdwvol=basinvdwvol+basinvdwvolp*dvol
else if (itype==10) then
	eleint=eleint+eleintp*dvol
	xint=xint+xintp*dvol
	yint=yint+yintp*dvol
	zint=zint+zintp*dvol
	xxint=xxint+xxintp*dvol
	yyint=yyint+yyintp*dvol
	zzint=zzint+zzintp*dvol
	xyint=xyint+xyintp*dvol
	yzint=yzint+yzintp*dvol
	xzint=xzint+xzintp*dvol
end if
!$OMP end CRITICAL
!$OMP END PARALLEL

10 continue
!---- Exact refinement with/without multi-level splitting of boundary grids
if (itype==2.or.itype==3) then
	if (itype==3) then !Calculate grid data of gradient of electron density used to linear interpolation to obtain the value at any point
		call gengradmat
	end if
	write(*,*) "Integrating grids at basin boundary..." 
	nrk4lim=100
	nrk4gradswitch=40
	hsizeinit=0.25D0
	ifinish=0
	!$OMP PARALLEL private(ix,iy,iz,iatt,rnowx,rnowy,rnowz,rx,ry,rz,tmpval,tmpval2,tmpval3,intvalp,basinvolp,basinvdwvolp,dist,tmps,iter,switchwei,&
	!$OMP rnowxtmp,rnowytmp,rnowztmp,orgxref,orgyref,orgzref,dxref,dyref,dzref,ixref,iyref,izref,nrefine,ndiv,&
	!$OMP k1,k2,k3,k4,dens,denshold,grad,hess,iattref,xtmp,ytmp,ztmp,irk4,hsize,ixtest,iytest,iztest,tmpdist) shared(intval,basinvol,basinvdwvol,ifinish) NUM_THREADS(nthreads)
	intvalp=0D0
	basinvolp=0D0
	basinvdwvolp=0D0
	!$OMP do schedule(DYNAMIC)
	do iz=2,nz-1
		rnowz=zarr(iz)
		do iy=2,ny-1
			rnowy=yarr(iy)
			do ix=2,nx-1
				rnowx=xarr(ix)
				if (.not.interbasgrid(ix,iy,iz)) cycle
! 				if (cubmat(ix,iy,iz)>0.001D0) then
! 					nrefine=2 !3 is the best
! 				else if (cubmat(ix,iy,iz)>0.001D0) then !0.0001 is the best
! 					nrefine=1
! 				else
! 					nrefine=1
! 				end if
 				nrefine=1
				ndiv=nrefine**3
				orgxref=rnowx-dx/2 !Take corner position as original point of microcycle
				orgyref=rnowy-dy/2
				orgzref=rnowz-dz/2
				dxref=dx/nrefine
				dyref=dy/nrefine
				dzref=dz/nrefine
				do ixref=1,nrefine
					do iyref=1,nrefine
						do izref=1,nrefine
							rnowxtmp=orgxref+(ixref-0.5D0)*dxref !Coordinate of current refined grid
							rnowytmp=orgyref+(iyref-0.5D0)*dyref
							rnowztmp=orgzref+(izref-0.5D0)*dzref
							if (cubmat(ix,iy,iz)<=0.001D0) then !Only refine the boundary inside vdW surface
								iattref=gridbas(ix,iy,iz)
							else
								xtmp=rnowxtmp !This point will continuously move in the iteration
								ytmp=rnowytmp
								ztmp=rnowztmp
								hsize=hsizeinit
								densold=0D0
								!** Tracing steepest ascent trajectory using 4-order Runge-Kutta (RK4)
		cycrk4:					do irk4=1,nrk4lim
									!For full accuracy refinement, or the first step, or when interpolation gradient works worse,&
									!namely has not converge until nrk4gradswitch, use exactly evaluated gradient
									if (itype==2.or.irk4==1.or.irk4==2.or.irk4>nrk4gradswitch) then 
										if (itype==3.and.irk4==nrk4gradswitch+1) then !Interpolated gradient doesn't work well, switch to full accuracy, reset the coordinate
											xtmp=rnowxtmp
											ytmp=rnowytmp
											ztmp=rnowztmp
											hsize=hsizeinit
										end if
										call calchessmat_dens(1,xtmp,ytmp,ztmp,dens,grad,hess) !Only value and gradient
										if (dens<densold-1D-10) then
											hsize=hsize*0.75D0 !Reduce step size if density decrease
										else if (dens>densold+1D-10) then
											hsize=hsizeinit !Recover to initial step size
										end if
										denshold=dens
										k1=grad/dsqrt(sum(grad**2))
										call calchessmat_dens(1,xtmp+hsize/2*k1(1),ytmp+hsize/2*k1(2),ztmp+hsize/2*k1(3),dens,grad,hess) !Only value and gradient
										k2=grad/dsqrt(sum(grad**2))
										call calchessmat_dens(1,xtmp+hsize/2*k2(1),ytmp+hsize/2*k2(2),ztmp+hsize/2*k2(3),dens,grad,hess) !Only value and gradient
										k3=grad/dsqrt(sum(grad**2))
										call calchessmat_dens(1,xtmp+hsize*k3(1),ytmp+hsize*k3(2),ztmp+hsize*k3(3),dens,grad,hess) !Only value and gradient
										k4=grad/dsqrt(sum(grad**2))
									else !Using the gradients evaluated by trilinear interpolation from pre-calculated grid data to save computational time
										call linintp3dvec(xtmp,ytmp,ztmp,grad) !Only value and gradient
! 										if (dens<densold-1D-10) then
! 											hsize=hsize*0.75D0
! 										else if (dens>densold+1D-10) then
! 											hsize=hsizeinit !Recover to initial step size
! 										end if
! 										denshold=dens
										k1=grad/dsqrt(sum(grad**2))
										call linintp3dvec(xtmp+hsize/2*k1(1),ytmp+hsize/2*k1(2),ztmp+hsize/2*k1(3),grad)
										k2=grad/dsqrt(sum(grad**2))
										call linintp3dvec(xtmp+hsize/2*k2(1),ytmp+hsize/2*k2(2),ztmp+hsize/2*k2(3),grad)
										k3=grad/dsqrt(sum(grad**2))
										call linintp3dvec(xtmp+hsize*k3(1),ytmp+hsize*k3(2),ztmp+hsize*k3(3),grad)
										k4=grad/dsqrt(sum(grad**2))
									end if
									xtmp=xtmp+hsize/6*(k1(1)+2*k2(1)+2*k3(1)+k4(1)) !Update current coordinate
									ytmp=ytmp+hsize/6*(k1(2)+2*k2(2)+2*k3(2)+k4(2))
									ztmp=ztmp+hsize/6*(k1(3)+2*k2(3)+2*k3(3)+k4(3))
									!Check if current position has entered trust radius of an attractor
									do iatttmp=1,numrealatt
										dist=dsqrt( (xtmp-CPpos(1,iatttmp))**2+(ytmp-CPpos(2,iatttmp))**2+(ztmp-CPpos(3,iatttmp))**2 )
										if (dist<trustrad(iatttmp)) then
											iattref=iatttmp
											exit cycrk4
										end if
									end do
									!Check if the closest grid and its 26 neighbours have the same attribution, if yes, employ its attribution then exit
									do ixtest=2,nx-1
										tmpdist=abs(xtmp-xarr(ixtest))
										if (tmpdist<dx/2D0) exit
									end do
									do iytest=2,ny-1
										tmpdist=abs(ytmp-yarr(iytest))
										if (tmpdist<dy/2D0) exit
									end do
									do iztest=2,nz-1
										tmpdist=abs(ztmp-zarr(iztest))
										if (tmpdist<dz/2D0) exit
									end do
									iattref=gridbas(ixtest,iytest,iztest)
									do imove=1,26
										if ( gridbas(ixtest+vec26x(imove),iytest+vec26y(imove),iztest+vec26z(imove))/=iattref ) exit
									end do
									if (imove==27) exit !Successfully passed neighbour test
								end do cycrk4
								if (irk4==nrk4lim+1) then !Didn't enter trust radius or didn't approach a grid who and whose neighbour have the same attribution
! 									write(*,*) "Warning: Exceeded the step limit of steepest ascent process!" !may frighten users, so comment
									iattref=gridbas(ix,iy,iz) !Use its original attribution
								end if
! 								write(*,*) irk4
							end if
							gridbas(ix,iy,iz)=iattref !Update attribution of boundary grids
							!Calculate switching function at current grid
							rx=rnowxtmp-CPpos(1,iattref) !The relative distance between current point to corresponding attractor
							ry=rnowytmp-CPpos(2,iattref)
							rz=rnowztmp-CPpos(3,iattref)
							dist=dsqrt(rx*rx+ry*ry+rz*rz)
							tmps=dist-trustrad(iattref)
							if (tmps>1) then
								switchwei=0
							else if (tmps<-1) then
								switchwei=1
							else
								do iter=1,nbeckeiter
									tmps=1.5D0*(tmps)-0.5D0*(tmps)**3
								end do
								switchwei=0.5D0*(1-tmps)
							end if
							switchwei=1-switchwei
							basinvolp(iattref)=basinvolp(iattref)+1D0/ndiv !Calculate boundary basin volume
							if (cubmat(ix,iy,iz)>0.001D0) basinvdwvolp(iattref)=basinvdwvolp(iattref)+1D0/ndiv
							if (ispecial==2.or.ifuncint==-1) then
								continue !Don't calculate function value, but only update attribution of boundary grids at this stage
							else if (ispecial==0) then
								tmpval=calcfuncall(ifuncint,rnowxtmp,rnowytmp,rnowztmp)
								intvalp(iattref,1)=intvalp(iattref,1)+tmpval*switchwei/ndiv
							else if (ispecial==1) then
								tmpval=infoentro(2,rnowxtmp,rnowytmp,rnowztmp) !Shannon entropy density, see JCP,126,191107 for example
								tmpval2=Fisherinfo(1,rnowxtmp,rnowytmp,rnowztmp) !Fisher information density, see JCP,126,191107 for example
								tmpval3=weizsacker(rnowxtmp,rnowytmp,rnowztmp) !Steric energy
								tmpval4=fdens(rnowxtmp,rnowytmp,rnowztmp) !Electron density
								intvalp(iattref,1)=intvalp(iattref,1)+tmpval*switchwei/ndiv
								intvalp(iattref,2)=intvalp(iattref,2)+tmpval2*switchwei/ndiv
								intvalp(iattref,3)=intvalp(iattref,3)+tmpval3*switchwei/ndiv
								intvalp(iattref,4)=intvalp(iattref,4)+tmpval4*switchwei/ndiv
							end if
						end do !End refine grid
					end do
				end do
				
			end do !End cycle ix grid
		end do
		ifinish=ifinish+1
		call showprog(ifinish,nz-2)
	end do
	!$OMP end do
	!$OMP CRITICAL
	intval=intval+intvalp*dvol
	basinvol=basinvol+basinvolp*dvol
	basinvdwvol=basinvdwvol+basinvdwvolp*dvol
	!$OMP end CRITICAL
	!$OMP END PARALLEL
	call detectinterbasgrd(6)
	write(*,*) "Basin boundary has been updated"
	numinterbas=count(interbasgrid)
	write(*,"(' The number of interbasin grids:',i12)") numinterbas	
end if

!Below are special modules for the cases when density of atom in free-state are involved
if (ispecial==2) then !Shubin's 2nd project, integrate relative Shannon and relative Fisher information
	allocate(rhogrid(nx,ny,nz),rhogradgrid(3,nx,ny,nz))
	ifinish=0
	write(*,*)
	write(*,*) "Calculating electron density and its gradient for actual system at each grid"
	!$OMP PARALLEL DO SHARED(rhogrid,rhogradgrid,ifinish) PRIVATE(ix,iy,iz,ptx,pty,ptz) schedule(dynamic) NUM_THREADS(nthreads)
	do iz=2,nz-1
		ptz=zarr(iz)
		do iy=2,ny-1
			pty=yarr(iy)
			do ix=2,nx-1
				ptx=xarr(ix)
				call calchessmat_dens(1,ptx,pty,ptz,rhogrid(ix,iy,iz),rhogradgrid(:,ix,iy,iz),hess)
			end do
		end do
		!$OMP CRITICAL
		ifinish=ifinish+1
		write(*,"(' Finished',i6,' /',i6)") ifinish,nz-1
		!$OMP end CRITICAL
	end do
	!$OMP end PARALLEL DO
	write(*,*)
	write(*,*) "Calculating electron density and its gradient for free-state atom at each grid"
	do iatt=1,numrealatt !Cycle each attractors
		write(*,"(' Processing ',a)") trim(custommapname(att2atm(iatt)))
		call dealloall
		call readwfn(custommapname(att2atm(iatt)),1)
		!$OMP PARALLEL private(intvalp,ix,iy,iz,ptx,pty,ptz,rx,ry,rz,dist,tmps,switchwei,&
		!$OMP rho,rhograd2,prodens,prodensgrad,prodensgrad2,tmpval1,tmpval2,tmpval3,tmpx,tmpy,tmpz) shared(intval) NUM_THREADS(nthreads)
		intvalp=0D0
		!$OMP do schedule(DYNAMIC)
		do iz=2,nz-1
			ptz=zarr(iz)
			do iy=2,ny-1
				pty=yarr(iy)
				do ix=2,nx-1
					ptx=xarr(ix)
					if (gridbas(ix,iy,iz)==iatt) then
						!Calculate switching function at current grid
						rx=ptx-CPpos(1,iatt) !The relative distance between current point to corresponding attractor
						ry=pty-CPpos(2,iatt)
						rz=ptz-CPpos(3,iatt)
						dist=dsqrt(rx*rx+ry*ry+rz*rz)
						tmps=dist-trustrad(iatt)
						if (tmps>1) then
							switchwei=0
						else if (tmps<-1) then
							switchwei=1
						else
							do iter=1,nbeckeiter
								tmps=1.5D0*(tmps)-0.5D0*(tmps)**3
							end do
							switchwei=0.5D0*(1-tmps)
						end if
						switchwei=1-switchwei
						rho=rhogrid(ix,iy,iz)
						rhograd2=sum(rhogradgrid(1:3,ix,iy,iz)**2)
						call calchessmat_dens(1,ptx,pty,ptz,prodens,prodensgrad(1:3),hess)
						prodensgrad2=sum(prodensgrad(1:3)**2)
						!Relative Shannon entropy
						tmpval1=rho*log(rho/prodens)
						intvalp(iatt,1)=intvalp(iatt,1)+tmpval1*switchwei
						!Relative Fisher information (old and incorrect form)
						tmpval2=rhograd2/rho-prodensgrad2/prodens
						intvalp(iatt,2)=intvalp(iatt,2)+tmpval2*switchwei
						!Relative Fisher information (new and correct form)
						tmpx=rhogradgrid(1,ix,iy,iz)/rho-prodensgrad(1)/prodens
						tmpy=rhogradgrid(2,ix,iy,iz)/rho-prodensgrad(2)/prodens
						tmpz=rhogradgrid(3,ix,iy,iz)/rho-prodensgrad(3)/prodens
						tmpval3=(tmpx**2+tmpy**2+tmpz**2)*rho
						intvalp(iatt,3)=intvalp(iatt,3)+tmpval3*switchwei
					end if
				end do
			end do
		end do
		!$OMP end do
		!$OMP CRITICAL
		intval=intval+intvalp*dvol
		!$OMP end CRITICAL
		!$OMP END PARALLEL
	end do
	deallocate(rhogrid,rhogradgrid)
	call dealloall
	write(*,"(' Reloading ',a)") trim(firstfilename)
	call readinfile(firstfilename,1) !Retrieve to first loaded file(whole molecule)
else if (ifuncint==-1) then !Deformation density
	allocate(prorhogrid(nx,ny,nz))
	write(*,*)
	write(*,*) "Calculating promolecular density at each grid"
	prorhogrid=0D0
	do iatm=1,ncenter_org !Cycle each atom
		write(*,"(' Processing atom',i6,a,'...')") iatm,a_org(iatm)%name
		call dealloall
		call readwfn(custommapname(iatm),1)
		!$OMP PARALLEL DO SHARED(prorhogrid) PRIVATE(ix,iy,iz) schedule(dynamic) NUM_THREADS(nthreads)
		do iz=2,nz-1
			do iy=2,ny-1
				do ix=2,nx-1
					prorhogrid(ix,iy,iz)=prorhogrid(ix,iy,iz)+fdens(xarr(ix),yarr(iy),zarr(iz))
				end do
			end do
		end do
		!$OMP end PARALLEL DO
	end do
	do iatt=1,numrealatt !Cycle each attractors
		do iz=2,nz-1
			do iy=2,ny-1
				do ix=2,nx-1
					if (gridbas(ix,iy,iz)==iatt) then
						!Calculate switching function at current grid
						rx=xarr(ix)-CPpos(1,iatt) !The relative distance between current point to corresponding attractor
						ry=yarr(iy)-CPpos(2,iatt)
						rz=zarr(iz)-CPpos(3,iatt)
						dist=dsqrt(rx*rx+ry*ry+rz*rz)
						tmps=dist-trustrad(iatt)
						if (tmps>1) then
							switchwei=0
						else if (tmps<-1) then
							switchwei=1
						else
							do iter=1,nbeckeiter
								tmps=1.5D0*(tmps)-0.5D0*(tmps)**3
							end do
							switchwei=0.5D0*(1-tmps)
						end if
						switchwei=1-switchwei
						defdens=cubmat(ix,iy,iz)-prorhogrid(ix,iy,iz)
						intval(iatt,1)=intval(iatt,1)+defdens*switchwei*dvol
					end if
				end do
			end do
		end do
	end do
	deallocate(prorhogrid)
	call dealloall
	call readinfile(firstfilename,1) !Retrieve to first loaded file(whole molecule) to calc real rho again
end if

!!----------- Output SUMMARY
if (itype==1.or.itype==2.or.itype==3) then !Integrate specific real space function(s)
	write(*,*)
	if (ifuncint==1) then
		write(*,*) "Total result:"
		write(*,*) "  #Basin        Integral(a.u.)      Vol(Bohr^3)    Vol(rho>0.001)"
		do iatt=1,numrealatt
			write(*,"(i8,f22.10,2f16.3)") iatt,intval(iatt,1),basinvol(iatt),basinvdwvol(iatt)
		end do
		write(*,"(' Sum of above integrals:',f20.8)") sum(intval(1:numrealatt,1))
		write(*,"(' Sum of basin volumes (rho>0.001):',f12.3,' Bohr^3')") sum(basinvdwvol(1:numrealatt))
		if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==0)) write(*,"(' Integral of unassigned grids:',f20.8)") intval(0,1)
		if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==-1)) write(*,"(' Integral of the grids travelled to box boundary:',f20.8)") intval(-1,1)
		write(*,*)
		rnormfac=sum(intval(1:numrealatt,1))/(nelec+nEDFelec) !The electrons represented by EDF must be taken into account!
		write(*,"(' Normalization factor of the integral of electron density is',f12.6)") rnormfac
		write(*,*) "The atomic charges after normalization and atomic volumes:"
		do iatm=0,ncenter
			do iatt=1,numrealatt
				if (att2atm(iatt)==iatm.and.iatm==0) then
					write(*,"(i7,' (NNA)   Charge:',f12.6,'     Volume:',f10.3,' Bohr^3')") iatt,-intval(iatt,1)/rnormfac,basinvdwvol(iatt)
				else if (att2atm(iatt)==iatm.and.iatm/=0) then
				    if (nEDFelec==0) then !Normal case, all electron basis or using pseudopotential but not accompanied by EDF
					    write(*,"(i7,' (',a,')    Charge:',f12.6,'     Volume:',f10.3,' Bohr^3')") iatm,a(iatm)%name,a(iatm)%charge-intval(iatt,1)/rnormfac,basinvdwvol(iatt)
					else !EDF is used, so using a(iatm)%index instead of a(iatm)%charge
					    write(*,"(i7,' (',a,')    Charge:',f12.6,'     Volume:',f10.3,' Bohr^3')") iatm,a(iatm)%name,a(iatm)%index-intval(iatt,1)/rnormfac,basinvdwvol(iatt)
					end if
				end if
			end do
		end do
	else
		if (ispecial==0) then
			write(*,*) "Total result:"
			write(*,*) "    Atom       Basin       Integral(a.u.)   Vol(Bohr^3)   Vol(rho>0.001)"
			do iatm=0,ncenter
				do iatt=1,numrealatt
					if (att2atm(iatt)==iatm.and.iatm==0) then
						write(*,"('      NNA   ',i8,f20.8,2f14.3)") iatt,intval(iatt,1),basinvol(iatt),basinvdwvol(iatt)
					else if (att2atm(iatt)==iatm.and.iatm/=0) then
						write(*,"(i7,' (',a,')',i8,f20.8,2f14.3)") iatm,a(iatm)%name,iatt,intval(iatt,1),basinvol(iatt),basinvdwvol(iatt)
					end if
				end do
			end do
			write(*,"(' Sum of above integrals:',f23.8)") sum(intval(1:numrealatt,1))
			write(*,"(' Sum of basin volumes (rho>0.001):',f12.3,' Bohr^3')") sum(basinvdwvol(1:numrealatt))
			if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==0)) write(*,"(' Integral of unassigned grids:',f20.8)") intval(0,1)
			if (any(gridbas(2:nx-1,2:ny-1,2:nz-1)==-1)) write(*,"(' Integral of the grids travelled to box boundary:',f20.8)") intval(-1,1)
		else if (ispecial==1) then
			write(*,*) "Total result:"
			write(*,*) "    Atom       Basin     Shannon       Fisher      Steric ene  Vol(rho>0.001)"
			do iatm=0,ncenter
				do iatt=1,numrealatt
					if (att2atm(iatt)==iatm.and.iatm==0) then
						write(*,"('      NNA   ',i8,3f14.7,f13.5)") iatt,intval(iatt,1),intval(iatt,1:3),basinvdwvol(iatt)
					else if (att2atm(iatt)==iatm.and.iatm/=0) then
						write(*,"(i7,' (',a,')',i8,3f14.7,f13.5)") iatm,a(iatm)%name,iatt,intval(iatt,1:3),basinvdwvol(iatt)
					end if
				end do
			end do
			write(*,"(' Total Shannon entropy:   ',f23.8)") sum(intval(1:numrealatt,1))
			write(*,"(' Total Fisher information:',f23.8)") sum(intval(1:numrealatt,2))
			write(*,"(' Total steric energy:     ',f23.8)") sum(intval(1:numrealatt,3))
			write(*,"(' Sum of basin volumes (rho>0.001):',f12.3,' Bohr^3')") sum(basinvdwvol(1:numrealatt))
			rnormfac=sum(intval(1:numrealatt,4))/(nelec+nEDFelec) !The electrons represented by EDF must be taken into account!
			write(*,"(/,' Normalization factor of the integral of electron density is',f12.6)") rnormfac
			write(*,*) "The atomic charges after normalization and atomic volumes:"
			do iatm=0,ncenter
				do iatt=1,numrealatt
					if (att2atm(iatt)==iatm.and.iatm==0) then
						write(*,"(i7,' (NNA)   Charge:',f12.6,'     Volume:',f10.3,' Bohr^3')") iatt,-intval(iatt,1)/rnormfac,basinvdwvol(iatt)
					else if (att2atm(iatt)==iatm.and.iatm/=0) then
						if (nEDFelec==0) then !Normal case, all electron basis or using pseudopotential but not accompanied by EDF
							write(*,"(i7,' (',a,')    Charge:',f12.6,'     Volume:',f10.3,' Bohr^3')") iatm,a(iatm)%name,a(iatm)%charge-intval(iatt,4)/rnormfac,basinvdwvol(iatt)
						else !EDF is used, so using a(iatm)%index instead of a(iatm)%charge
							write(*,"(i7,' (',a,')    Charge:',f12.6,'     Volume:',f10.3,' Bohr^3')") iatm,a(iatm)%name,a(iatm)%index-intval(iatt,4)/rnormfac,basinvdwvol(iatt)
						end if
					end if
				end do
			end do
		else if (ispecial==2) then
			write(*,*) "Total result:"
			write(*,*) "    Atom      Basin        Rel.Shannon     Rel.Fisher(old)     Rel.Fisher(new)"
			do iatm=1,ncenter
				do iatt=1,numrealatt
					if (att2atm(iatt)==iatm) write(*,"(i7,' (',a,')',i7,3f20.8)") iatm,a(iatm)%name,iatt,intval(iatt,1:3)
				end do
			end do
			write(*,"(' Sum of relat_Shannon:     ',f23.8)") sum(intval(1:numrealatt,1))
			write(*,"(' Sum of relat_Fisher(old): ',f23.8)") sum(intval(1:numrealatt,2))
			write(*,"(' Sum of relat_Fisher(new): ',f23.8)") sum(intval(1:numrealatt,3))
		end if
	end if
	write(*,*)
else if (itype==10) then !Electric multipole moment
	ioutid=6
101	eleinttot=0D0
	dipelextot=0D0
	dipeleytot=0D0
	dipeleztot=0D0
	if (ioutid==6) write(*,*)
	write(ioutid,*) "Note: All units shown below are in a.u.!"
	write(ioutid,*)
	do iatm=0,ncenter
		do iatt=1,numrealatt
			if (att2atm(iatt)/=iatm) cycle
			if (iatm==0) then
				write(ioutid,"(' Result of NNA',i6)") iatt
			else
				write(ioutid,"(' Result of atom',i6,' (',a,')')") iatm,a(iatm)%name
			end if
			write(ioutid,"(' Basin electric monopole moment:',f12.6)") -eleint(iatt)
			write(ioutid,"(' Basin electric dipole moment:')") 
			write(ioutid,"('  X=',f12.6,'  Y=',f12.6,'  Z=',f12.6,'  Magnitude=',f12.6)") -xint(iatt),-yint(iatt),-zint(iatt),sqrt(xint(iatt)**2+yint(iatt)**2+zint(iatt)**2)
			eleinttot=eleinttot+eleint(iatt)
			dipelex=-eleint(iatt)*CPpos(1,iatt)+(-xint(iatt)) !Contribution to molecular total dipole moment
			dipeley=-eleint(iatt)*CPpos(2,iatt)+(-yint(iatt))
			dipelez=-eleint(iatt)*CPpos(3,iatt)+(-zint(iatt))
			dipelextot=dipelextot+dipelex
			dipeleytot=dipeleytot+dipeley
			dipeleztot=dipeleztot+dipelez
			write(ioutid,"(' Basin electron contribution to molecular dipole moment:')")
			write(ioutid,"('  X=',f12.6,'  Y=',f12.6,'  Z=',f12.6,'  Magnitude=',f12.6)") dipelex,dipeley,dipelez,sqrt(dipelex**2+dipeley**2+dipelez**2)
			write(ioutid,"(' Basin electric quadrupole moment (Cartesian form):')")
			rrint=xxint(iatt)+yyint(iatt)+zzint(iatt)
			QXX=-(3*xxint(iatt)-rrint)/2
			QYY=-(3*yyint(iatt)-rrint)/2
			QZZ=-(3*zzint(iatt)-rrint)/2
			write(ioutid,"(' QXX=',f12.6,'  QXY=',f12.6,'  QXZ=',f12.6)") QXX,-(3*xyint(iatt))/2,-(3*xzint(iatt))/2
			write(ioutid,"(' QYX=',f12.6,'  QYY=',f12.6,'  QYZ=',f12.6)") -(3*xyint(iatt))/2,QYY,-(3*yzint(iatt))/2
			write(ioutid,"(' QZX=',f12.6,'  QZY=',f12.6,'  QZZ=',f12.6)") -(3*xzint(iatt))/2,-(3*yzint(iatt))/2,QZZ
			write(ioutid,"( ' The magnitude of electric quadrupole moment (Cartesian form):',f12.6)") dsqrt(2D0/3D0*(QXX**2+QYY**2+QZZ**2))
			R20=-(3*zzint(iatt)-rrint)/2D0 !Notice that the negative sign, because electrons carry negative charge
			R2n1=-dsqrt(3D0)*yzint(iatt)
			R2p1=-dsqrt(3D0)*xzint(iatt)
			R2n2=-dsqrt(3D0)*xyint(iatt)
			R2p2=-dsqrt(3D0)/2D0*(xxint(iatt)-yyint(iatt))
			write(ioutid,"(' Electric quadrupole moments (Spherical harmonic form):')")
			write(ioutid,"(' Q_2,0 =',f11.6,'   Q_2,-1=',f11.6,'   Q_2,1=',f11.6)") R20,R2n1,R2p1
			write(ioutid,"(' Q_2,-2=',f11.6,'   Q_2,2 =',f11.6)") R2n2,R2p2
			write(ioutid,"( ' Magnitude: |Q_2|=',f12.6)") dsqrt(R20**2+R2n1**2+R2p1**2+R2n2**2+R2p2**2)
			write(ioutid,*)
		end do
	end do
	!Output overall electric properties
	dipnucx=sum(a(:)%x*a(:)%charge)
	dipnucy=sum(a(:)%y*a(:)%charge)
	dipnucz=sum(a(:)%z*a(:)%charge)
	write(ioutid,"( ' Molecular net charge:',f12.6)") sum(a%charge)-eleinttot
	write(ioutid,"( ' Nuclear contribution to molecular dipole moment:')") 
	write(ioutid,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") dipnucx,dipnucy,dipnucz,sqrt(dipnucx**2+dipnucy**2+dipnucz**2)
	write(ioutid,"( ' Electron contribution to molecular dipole moment:')") 
	write(ioutid,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") dipelextot,dipeleytot,dipeleztot,sqrt(dipelextot**2+dipeleytot**2+dipeleztot**2)
	dipmolx=dipnucx+dipelextot
	dipmoly=dipnucy+dipeleytot
	dipmolz=dipnucz+dipeleztot
	write(ioutid,"( ' Molecular dipole moment:')")
	write(ioutid,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") dipmolx,dipmoly,dipmolz,sqrt(dipmolx**2+dipmoly**2+dipmolz**2)
	if (ioutid==6) write(*,*)
	if (ioutid==10) then
		close(10)
		write(*,*) "Done!"
		return
	end if
end if

CALL CPU_TIME(time_end)
call walltime(walltime2)
write(*,"(' Integrating basins took up CPU time',f12.2,'s, wall clock time',i10,'s')") time_end-time_begin,walltime2-walltime1

if (itype==10) then
	write(*,*)
	write(*,*) "If also outputting the result to multipol.txt in current folder? (y/n)"
	read(*,*) selectyn
	if (selectyn=='y'.or.selectyn=='Y') then
		ioutid=10
		open(10,file="multipol.txt",status="replace")
		goto 101
	end if
end if
end subroutine





!!------ Generate grid data of gradient of electron density, used to refine basin boundary
subroutine gengradmat
use defvar
use util
use function
implicit real*8 (a-h,o-z)
real*8 wfnval(nmo),wfnderv(3,nmo),gradrho(3),EDFgrad(3),sumgrad2
if (allocated(cubmatvec)) then
	if (size(cubmatvec,2)==nx.and.size(cubmatvec,3)==ny.and.size(cubmatvec,4)==nz) return !Already generated
	deallocate(cubmatvec)
end if
allocate(cubmatvec(3,nx,ny,nz))
ifinish=0
write(*,*) "Generating grid data of gradient of electron density..."
!$OMP PARALLEL DO SHARED(cubmatvec,ifinish) PRIVATE(ix,iy,iz,tmpx,tmpy,tmpz,wfnval,wfnderv,gradrho,imo) schedule(dynamic) NUM_THREADS(nthreads)
do iz=1,nz
	tmpz=orgz+(iz-1)*dz
	do iy=1,ny
		tmpy=orgy+(iy-1)*dy
		do ix=1,nx
			tmpx=orgx+(ix-1)*dx
			call orbderv(2,1,nmo,tmpx,tmpy,tmpz,wfnval,wfnderv)
			gradrho=0D0
			do imo=1,nmo
				gradrho(:)=gradrho(:)+MOocc(imo)*wfnval(imo)*wfnderv(:,imo)
			end do
			cubmatvec(:,ix,iy,iz)=2*gradrho(:)
		end do
	end do
    ifinish=ifinish+1
    call showprog(ifinish,nz)
end do
!$OMP END PARALLEL DO
end subroutine



!!------- Calculate high ELF localization domain populations and volumes (HELP and HELV), ChemPhysChem, 14, 3714 (2013)
subroutine HELP_HELV
use defvar
use function
use basinintmod
use util
implicit real*8(a-h,o-z)
real*8 :: ELFthres=0.5D0,rhothres=0.001D0
real*8 HELP,HELV,HELP_priv,HELV_priv
character c2000tmp*2000
integer tmparr(numrealatt)

do while(.true.)
    write(*,*)
    write(*,*) "    -------------------- Calculation of HELP and HELV --------------------"
    write(*,*) "-1 Return"
    write(*,*) " 0 Select basins and calculate their HELP and HELV"
    write(*,"(a,f6.3)") "  1 Define threshold of ELF, current: >",ELFthres
    write(*,"(a,f8.5,' a.u.')") "  2 Define threshold of electron density, current: >",rhothres
    read(*,*) isel
    if (isel==-1) then
        return
    else if (isel==1) then
        write(*,*) "Input threshold of ELF, e.g. 0.5"
        read(*,*) ELFthres
    else if (isel==2) then
        write(*,*) "Input threshold of electron density in a.u., e.g. 0.001"
        read(*,*) rhothres
    else if (isel==0) then
        write(*,*) "Input indices of the ELF basins for which HELP and HELV will be calculated"
        write(*,*) "e.g. 2,3,7-10,14"
        read(*,"(a)") c2000tmp
        call str2arr(c2000tmp,nselbasin,tmparr)

        write(*,*) "Integrating, please wait patiently..."
        do idx=1,nselbasin
            ibasin=tmparr(idx)
            HELP=0
            HELV=0
            basinpop=0
            basinvol=0
            !$OMP PARALLEL private(ix,iy,iz,irealatt,rnowx,rnowy,rnowz,HELP_priv,HELV_priv,basinpop_priv,basinvol_priv,rhotmp,ELFtmp) &
            !$OMP shared(HELP,HELV,basinpop,basinvol) NUM_THREADS(nthreads)
            HELP_priv=0
            HELV_priv=0
            basinpop_priv=0
            basinvol_priv=0
            !$OMP do schedule(DYNAMIC)
            do iz=2,nz-1
	            do iy=2,ny-1
		            do ix=2,nx-1
			            rnowx=orgx+(ix-1)*dx
			            rnowy=orgy+(iy-1)*dy
			            rnowz=orgz+(iz-1)*dz
			            irealatt=gridbas(ix,iy,iz)
                        rhotmp=fdens(rnowx,rnowy,rnowz)
                        ELFtmp=ELF_LOL(rnowx,rnowy,rnowz,"ELF")
                        if (irealatt==ibasin) then
                            if (rhotmp>rhothres.and.ELFtmp>ELFthres) then
                                HELP_priv=HELP_priv+rhotmp
                                HELV_priv=HELV_priv+1
                            end if
                            basinpop_priv=basinpop_priv+rhotmp
                            basinvol_priv=basinvol_priv+1
                        end if
		            end do
	            end do
            end do
            !$OMP end do
            !$OMP CRITICAL
                HELP=HELP+HELP_priv
                HELV=HELV+HELV_priv
                basinpop=basinpop+basinpop_priv
                basinvol=basinvol+basinvol_priv
            !$OMP end CRITICAL
            !$OMP END PARALLEL
            dvol=dx*dy*dz
            HELP=HELP*dvol
            HELV=HELV*dvol
            basinpop=basinpop*dvol
            basinvol=basinvol*dvol
            write(*,*)
            if (nselbasin>1) write(*,"(a,i6)") " ------ Information of basin",ibasin
            write(*,*) "Basin information: (constraints are not taken into account)"
            write(*,"(' Population:',f8.4,'   Volume:',f10.4,' Bohr^3')") basinpop,basinvol
            write(*,*)
            write(*,"(' High ELF localization domain population (HELP):',f10.4)") HELP
            write(*,"(' High ELF localization domain volume (HELV):    ',f10.4,' Bohr^3')") HELV
        end do
    end if
end do

end subroutine



!!------------ Determine correspondence between attractor index and atom index
!The passed in array is att2atm(numrealatt)
!Returned arry att2atm(i)=k means real attractor i corresponds to atom k. If k=0, that means attractor i should be NNA
subroutine calc_att2atm(att2atm)
use defvar
use basinintmod
implicit real*8 (a-h,o-z)
integer att2atm(numrealatt)
att2atm=0
crit=0.3D0 !If distance between attractor and a nucleus is smaller than 0.3 Bohr, then the attractor will belong to the atom
do iatt=1,numrealatt !Cycle each attractors
	do iatm=1,ncenter
        disttest=dsqrt( (realattxyz(iatt,1)-a(iatm)%x)**2+(realattxyz(iatt,2)-a(iatm)%y)**2+(realattxyz(iatt,3)-a(iatm)%z)**2 )
        if (disttest<0.3D0) then 
	        att2atm(iatt)=iatm
            exit
        end if
    end do
end do
end subroutine