!------------- Modify & Check wavefunction
subroutine modwfn
use defvar
use util
implicit real*8 (a-h,o-z)
character seltmpc*10,selectyn,c1000tmp*1000,c2000tmp*2000
real*8 eigval(nbasis),eigvec(nbasis,nbasis),tmpmat(nbasis,nbasis)
real*8,allocatable :: tmparr(:)
integer orbarr(nmo)
integer,allocatable :: exclfragatm(:),tmparrint(:),idxsel(:)
character(len=3) :: orbtype(0:2)=(/ "A+B"," A "," B " /)
character symstr*6

write(*,*) "Note: ""GTF"" in this function refers to Gaussian type primitive function"
call updatenelec !Sometimes the naelec/nbelec/nelec recorded in the input file is not in line with sum of occupation numbers, so update
do while(.true.)
	write(*,*)
	write(*,*) "            ============ Modify & Check wavefunction ============ "
	write(*,"(' GTF:',i6,', Orb:',i6,', Atom:',i5,', A/B/tot ele:',3f10.3)") nprims,nmo,ncenter,naelec,nbelec,nelec
	if (ifragcontri/=1) write(*,*) "-4 Exclude contribution of some atoms to real space functions"
	if (ifragcontri/=1) write(*,*) "-3 Only retain contribution of some atoms to real space functions"
	write(*,*) "-1 Return"
	write(*,*) "0 Save the modified wavefunction to a new .wfn file"
	if (allocated(CObasa)) then
		write(*,*) "1 List all GTFs      2 List all basis function"
	else
		write(*,*) "1 List all GTFs"
	end if
	write(*,*) "3 List all orbitals                4 Print detail information of an orbital"
	if (allocated(CObasa)) write(*,*) "5 Print coefficient matrix in basis function"
	if (allocated(CObasa)) write(*,*) "6 Print density matrix in basis function"
	if (allocated(CObasa)) write(*,*) "7 Print various kinds of integral matrix between basis functions"
	write(*,*) "11 Exchange some information of two GTFs"
	write(*,*) "21 Set center of a GTF             22 Set type of a GTF"
	write(*,*) "23 Set exponent of a GTF           24 Set coefficient of a GTF in an orbital"
	if (allocated(CObasa)) then
		write(*,*) "25 Set coefficients of GTFs/basis functions that satisfied certain conditions"
	else
		write(*,*) "25 Set coefficients of GTFs that satisfied certain conditions"
	end if
	write(*,*) "26 Set occupation of some orbitals  27 Set type of some orbitals"
	write(*,*) "28 Set energy of some orbitals      29 Exchange information of two orbitals"
    write(*,*) "30 Exchange energies and occupation numbers for all orbitals"
	write(*,*) "31 Translate the system             32 Translate and duplicate the system"
	write(*,*) "33 Rotate wavefunction, namely X->Y, Y->Z, Z->X"
	if (imodwfn==0) write(*,*) "34 Set occupation number of inner orbitals to zero" !If occupation has been modified, don't do this to complicate things
	if (allocated(MOsym)) write(*,*) "35 Keep or discard orbital contributions according to irreducible rep."
	write(*,*) "36 Invert phase of some orbitals"
	write(*,*) "37 Split spatial orbitals as alpha and beta spin orbitals"
	read(*,*) isel
	
	if (isel==-1) then
		if (allocated(CObasa).and.imodwfn==1) then
			write(*,*) "Updating density matrix..."
			call genP
			write(*,*) "Density matrix has been updated"
		end if
        if (imodwfn==1) if_initlibreta=0 !LIBRETA should then be re-initialized
		exit
		
	else if (isel==-3.or.isel==-4) then
		if (allocated(fragatm)) deallocate(fragatm) !fragatm has been defined previously by default, fragatm contains all atoms
		if (isel==-3) then
			! "fragatm" is convertion relationship from fragment to the whole,
			! e.g. fragatm(4) is the actual atom index corresponding the 4th atom in fragment list
			write(*,"(a)") " Input atomic indices to define the fragment, e.g. 1,3-6,8,10-11 means atoms 1,3,4,5,6,8,10,11 will constitute the fragment"
			read(*,"(a)") c2000tmp
			call str2arr(c2000tmp,nfragatm)
			allocate(fragatm(nfragatm))
			call str2arr(c2000tmp,nfragatm,fragatm)
			call sort(fragatm,"val")
		else if (isel==-4) then
			write(*,*) "Input indices of the atoms to be excluded"
			write(*,*) "e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will be excluded"
			read(*,"(a)") c2000tmp
			call str2arr(c2000tmp,nexclatm)
			nfragatm=ncenter-nexclatm
			allocate(fragatm(nfragatm),exclfragatm(nexclatm))
			call str2arr(c2000tmp,nexclatm,exclfragatm)
			j=0
			do i=1,ncenter
				if (all(exclfragatm/=i)) then
					j=j+1
					fragatm(j)=i
				end if
			end do
		end if
		j=0
		do i=1,nprims
			if (any(fragatm==b(i)%center)) then
				j=j+1      !Move function in the fragment to head of list
				CO(:,j)=CO(:,i)
				b(j)=b(i)
			end if
		end do
		ifragcontri=1 !Fragment has been defined by users
		write(*,"(' Done,',i8,' GTFs have been discarded,',i8,' GTFs reserved')") nprims-j,j
		nprims=j !Cut list at j, all functions after j seem non exist
		if (isel==-4) deallocate(exclfragatm)

		!Modification of wavefunction has finished, now reduce size of b, CO... to current nprims and nmo to avoid potential problems
		if (allocated(b)) then !Only for input file contains wavefunctions
			call resizebynmo(nmo,nprims) !Reduce size of CO, MOene, MOocc, MOtype
			allocate(b_tmp(nprims))
			b_tmp(:)=b(1:nprims)
			deallocate(b)
			allocate(b(nprims))
			b=b_tmp
			deallocate(b_tmp)
		end if
        
	else if (isel==0) then
		call outwfn("new.wfn",1,1,10)
		write(*,*) "Wavefunction has been outputted to new.wfn in current folder"
	
	else if (isel==1) then
		do i=1,nprims
			write(*,"(i6,' Center:',i5,'(',a2,')','   Type: ',a,'   Exponent:',E16.7)") i,b(i)%center,a(b(i)%center)%name,GTFtype2name(b(i)%type),b(i)%exp
		end do
	
	else if (isel==2) then
		do i=1,nbasis
            if (isphergau==1) then
			    write(*,"(' Basis:',i5,'   Shell:',i5,'   Center:',i5,'(',a2,')   Type:',a)")&
                i,basshell(i),bascen(i),a(bascen(i))%name,GTFtype2name(bastype(i))
            else !The primstart/end is constructed w.r.t. Cartesian basis functions
			    write(*,"(' Basis:',i5,'   Shell:',i5,'   Center:',i5,'(',a2,')   Type:',a,'  GTF:',i6,' to',i6)")&
                i,basshell(i),bascen(i),a(bascen(i))%name,GTFtype2name(bastype(i)),primstart(i),primend(i)
            end if
		end do
	
	else if (isel==3) then
		write(*,*) "Basic information of all orbitals:"
		symstr=" "
		naorb=count(MOtype==1)
		do i=1,nmo
			if (allocated(MOsym)) symstr='('//MOsym(i)//')'
			if (wfntype==0.or.wfntype==2.or.wfntype==3) then
				write(*,"(' Orb:',i6,' Ene(au/eV):',f13.6,f13.4,' Occ:',f9.6,' Type:',a,1x,a)") &
				i,MOene(i),MOene(i)*au2eV,MOocc(i),orbtype(MOtype(i)),symstr
			else
				if (MOtype(i)==1) then
					write(*,"(i6,9x,' E(au/eV):',f12.5,f13.4,' Occ:',f9.6,' Typ:',a,1x,a)") &
					i,MOene(i),MOene(i)*au2eV,MOocc(i),orbtype(MOtype(i)),symstr
				else
					write(*,"(i6,' (',i6,')',' E(au/eV):',f12.5,f13.4,' Occ:',f9.6,' Typ:',a,1x,a)") &
					i,i-naorb,MOene(i),MOene(i)*au2eV,MOocc(i),orbtype(MOtype(i)),symstr
				end if
			end if
		end do
		if (any(MOtype==2)) write(*,"(a)") " Note: For beta orbitals, &
		the index in the parenthese shown above is the index counted from the first beta orbital"
		
	else if (isel==4) then
		write(*,*) "Input the orbital index, e.g. 12"
		read(*,*) i
		if (i<1.or.i>nmo) then
			write(*,"(' Invalid orbital index, should within range of',i5,' and ',i5)") 1,nmo
		else
			write(*,"(' Occupation number is ',f12.7,'     Energy is',f12.6,' Hartree')") MOocc(i),MOene(i)
			if (MOtype(i)==0) write(*,*) "This is a closed-shell orbital"
			if (MOtype(i)==1) write(*,*) "This is an alpha orbital"
			if (MOtype(i)==2) write(*,*) "This is a beta orbital"
			write(*,*)
			do j=1,nprims
				write(*,"(' GTF:',i6,' Cen:',i5,'(',a2,')',' Type: ',a,' Coeff:',1PE16.8,' Exp: ',1PE13.7)") &
				j,b(j)%center,a(b(j)%center)%name,GTFtype2name(b(j)%type),CO(i,j),b(j)%exp
			end do
			write(*,"(a,/)") " Note: The ""coeff."" are expansion coefficients of orbitals with respect to GTFs, including normalization constant"
			if (allocated(b)) then
				do j=1,nbasis
					if (MOtype(i)==0.or.MOtype(i)==1) covalue=CObasa(j,i)
					if (MOtype(i)==2) covalue=CObasb(j,i-nbasis)
					write(*,"(' Basis func:',i6,'  Cen:',i5,'(',a2,')',' Shell:',i5,' Type: ',a,' Coeff:',f12.8)") &
					j,bascen(j),a(bascen(j))%name,basshell(j),GTFtype2name(bastype(j)),covalue
				end do
			end if
			write(*,"(a,/)") " Note: The ""coeff."" are expansion coefficients of orbitals with respect to basis functions, which are normalized functions"
		end if
	
	else if (isel==5) then
		write(*,*) "0 Return"
		write(*,*) "1 Print on screen"
		write(*,*) "2 Print to Cmat.txt in current folder"
		read(*,*) iseltmp
		if (iseltmp==1.or.iseltmp==2) then
			if (iseltmp==1) ides=6
			if (iseltmp==2) then
				ides=10
				open(ides,file="Cmat.txt",status="replace")
			end if
			write(ides,*) "Note: (i,j) element means coefficient of ith basis function in jth orbital"
			if (wfntype==0.or.wfntype==2.or.wfntype==3) then
				call showmatgau(CObasa,"Coefficient matrix",0,fileid=ides)
			else if (wfntype==1.or.wfntype==4) then
				call showmatgau(CObasa,"Alpha coefficient matrix",0,fileid=ides)
				call showmatgau(CObasb,"Beta coefficient matrix",0,fileid=ides)
			end if
			if (iseltmp==2) then
				write(*,*) "Done! The matrix has been outputted to Cmat.txt in current folder"
				close(ides)
			end if
		end if
	
	else if (isel==6) then
		write(*,*) "0 Return"
		write(*,*) "1 Print on screen"
		write(*,*) "2 Print to Pmat.txt in current folder"
		read(*,*) iseltmp
		if (iseltmp==1.or.iseltmp==2) then
			if (iseltmp==1) ides=6
			if (iseltmp==2) then
				ides=10
				open(ides,file="Pmat.txt",status="replace")
			end if
			call showmatgau(Ptot,"Total density matrix",1,fileid=ides)
			sumt=0
			do i=1,nbasis
				sumt=sumt+Ptot(i,i)
			end do
			write(ides,"(' The trace of the density matrix:',f12.6)") sumt
			write(ides,"(' The trace of the density matrix multiplied by overlap matrix:',f12.6)") sum(Ptot*Sbas)
			if (wfntype==1.or.wfntype==2.or.wfntype==4) then
				suma=0
				sumb=0
				do i=1,nbasis
					suma=suma+Palpha(i,i)
					sumb=sumb+Pbeta(i,i)
				end do
				write(ides,*)
				call showmatgau(Palpha-Pbeta,"Spin density matrix",1,fileid=ides)
				write(ides,*)
				call showmatgau(Palpha,"Alpha density matrix",1,fileid=ides)
				write(ides,*)
				call showmatgau(Pbeta,"Beta density matrix",1,fileid=ides)
				write(ides,*)
				write(ides,"(' The trace of the alpha and beta density matrix:',2f12.6)") suma,sumb
				write(ides,"(' The trace of the density matrix multiplied by overlap matrix:',/,' Alpha:',f12.6,'   Beta:',f12.6)") sum(Palpha*Sbas),sum(Pbeta*Sbas)
			end if
			if (iseltmp==2) then
				write(*,*) "Done! The matrix has been outputted to Pmat.txt in current folder"
				close(ides)
			end if
		end if
	
	!Print various kinds of integral matrix between basis functions
	else if (isel==7) then
		write(*,*) "Print which kind of integral matrix?"
		write(*,*) "1 Overlap integral"
		write(*,*) "2 Electric dipole moment integral"
		write(*,*) "3 Magnetic dipole moment integral"
		write(*,*) "4 Velocity integral"
		write(*,*) "5 Kinetic energy integral"
		write(*,*) "6 Electric quadrupole moment integral"
		write(*,*) "7 Electric octopole moment integral"
		write(*,*) "8 Electric hexadecapole moment integral"
		read(*,*) imattype
        iout=2
        if (imattype<=5) then !Amount of quadrupole and octopole is too large, and only be printed to file
		    write(*,*) "Select destination of outputting"
		    write(*,*) "1 Print on screen"
		    write(*,*) "2 Print to intmat.txt in current folder"
		    read(*,*) iout
        end if
		if (iout==1) then
            ides=6
		else
			ides=10
			open(ides,file="intmat.txt",status="replace")
		end if
		if (imattype==1) then
			call showmatgau(Sbas,"Overlap matrix",1,fileid=ides)
            tmpmat=Sbas
			call diagsymat(Sbas,eigvec,eigval,ierror)
			write(ides,*)
			write(ides,*) "Eigenvalues:"
			write(ides,"(6f12.8)") eigval
            Sbas=tmpmat
		else if (imattype==2) then
			if (.not.allocated(Dbas)) then
                write(*,*) "Calculating the matrix..."
                call genDbas_curr
            end if
			write(ides,*)
			call showmatgau(Dbas(1,:,:),"Electric dipole moment matrix (X component)",1,fileid=ides)
			write(ides,*)
			call showmatgau(Dbas(2,:,:),"Electric dipole moment matrix (Y component)",1,fileid=ides)
			write(ides,*)
			call showmatgau(Dbas(3,:,:),"Electric dipole moment matrix (Z component)",1,fileid=ides)
		else if (imattype==3) then
			if (.not.allocated(Magbas)) then
                write(*,*) "Calculating the matrix..."
                call genMagbas_curr
            end if
			write(ides,*)
			call showmatgau(Magbas(1,:,:),"Magnetic dipole moment matrix (X component)",0,fileid=ides)
			write(ides,*)
			call showmatgau(Magbas(2,:,:),"Magnetic dipole moment matrix (Y component)",0,fileid=ides)
			write(ides,*)
			call showmatgau(Magbas(3,:,:),"Magnetic dipole moment matrix (Z component)",0,fileid=ides)
		else if (imattype==4) then
			if (.not.allocated(Velbas)) then
                write(*,*) "Calculating the matrix..."
                call genVelbas_curr
            end if
			write(ides,*)
			call showmatgau(Velbas(1,:,:),"Velocity matrix (X component)",0,fileid=ides)
			write(ides,*)
			call showmatgau(Velbas(2,:,:),"Velocity matrix (Y component)",0,fileid=ides)
			write(ides,*)
			call showmatgau(Velbas(3,:,:),"Velocity matrix (Z component)",0,fileid=ides)
        else if (imattype==5) then
			if (.not.allocated(Tbas)) then
                write(*,*) "Calculating the matrix..."
                call genTbas_curr
            end if
			write(ides,*)
			call showmatgau(Tbas(:,:),"Kinetic energy matrix",1,fileid=ides)
        else if (imattype==6) then
			if (.not.allocated(Quadbas)) then
                write(*,*) "Calculating the matrix..."
                call genMultipolebas_curr
            end if
			              call showmatgau(Quadbas(1,:,:),"Quadrupole moment matrix (XX component)",1,fileid=ides)
			write(ides,*);call showmatgau(Quadbas(2,:,:),"Quadrupole moment matrix (YY component)",1,fileid=ides)
			write(ides,*);call showmatgau(Quadbas(3,:,:),"Quadrupole moment matrix (ZZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(Quadbas(4,:,:),"Quadrupole moment matrix (XY component)",1,fileid=ides)
			write(ides,*);call showmatgau(Quadbas(5,:,:),"Quadrupole moment matrix (YZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(Quadbas(6,:,:),"Quadrupole moment matrix (XZ component)",1,fileid=ides)
        else if (imattype==7) then
			if (.not.allocated(Octobas)) then
                write(*,*) "Calculating the matrix..."
                call genMultipolebas_curr
            end if
			              call showmatgau(Octobas(1,:,:),"Octopole moment matrix (XXX component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(2,:,:),"Octopole moment matrix (YYY component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(3,:,:),"Octopole moment matrix (ZZZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(4,:,:),"Octopole moment matrix (YZZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(5,:,:),"Octopole moment matrix (XZZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(6,:,:),"Octopole moment matrix (XXZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(7,:,:),"Octopole moment matrix (YYZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(8,:,:),"Octopole moment matrix (XXY component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(9,:,:),"Octopole moment matrix (XYY component)",1,fileid=ides)
			write(ides,*);call showmatgau(Octobas(10,:,:),"Octopole moment matrix (XYZ component)",1,fileid=ides)
        else if (imattype==8) then
			if (.not.allocated(hexdebas)) then
                write(*,*) "Calculating the matrix..."
                call genMultipolebas_curr
            end if
			              call showmatgau(hexdebas(1,:,:),"Hexadecapole moment matrix (XXXX component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(2,:,:),"Hexadecapole moment matrix (YYYY component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(3,:,:),"Hexadecapole moment matrix (ZZZZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(4,:,:),"Hexadecapole moment matrix (XXXY component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(5,:,:),"Hexadecapole moment matrix (XXXZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(6,:,:),"Hexadecapole moment matrix (YYYX component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(7,:,:),"Hexadecapole moment matrix (YYYZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(8,:,:),"Hexadecapole moment matrix (ZZZX component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(9,:,:),"Hexadecapole moment matrix (ZZZY component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(10,:,:),"Hexadecapole moment matrix (XXYY component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(11,:,:),"Hexadecapole moment matrix (XXZZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(12,:,:),"Hexadecapole moment matrix (YYZZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(13,:,:),"Hexadecapole moment matrix (XXYZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(14,:,:),"Hexadecapole moment matrix (YYXZ component)",1,fileid=ides)
			write(ides,*);call showmatgau(hexdebas(15,:,:),"Hexadecapole moment matrix (ZZXY component)",1,fileid=ides)
		end if
		if (iout==2) then
			write(*,*) "Done! The matrix has been outputted to intmat.txt in current folder"
			close(ides)
		end if
		
	else if (isel==11) then
		write(*,*) "Exchange information of which two GTFs? Input their indices, e.g. 18,21"
		read(*,*) i,j
		write(*,*) "Exchange which information for the two GTFs?"
		write(*,*) "1 Exchange all properties"
		write(*,*) "2 Exchange center"
		write(*,*) "3 Exchange function type"
		write(*,*) "4 Exchange exponent"
		write(*,*) "5 Exchange orbital expansion coefficient"
		read(*,*) iswapcontent
		if (iswapcontent==1) call swapGTF(i,j,"all")
		if (iswapcontent==2) call swapGTF(i,j,"cen")
		if (iswapcontent==3) call swapGTF(i,j,"typ")
		if (iswapcontent==4) call swapGTF(i,j,"exp")
		if (iswapcontent==5) call swapGTF(i,j,"MO ")
		write(*,*) "Exchanging finished!"
	
	else if (isel==21) then
		write(*,*) "Input index of the GTF, e.g. 13"
		read(*,*) i
		write(*,*) "Input the index of the center that you want the function centered at, e.g. 5"
		read(*,*) j
		if (j<=ncenter.and.j>0) then
			b(i)%center=j
            write(*,*) "Done!"
		else
			write(*,"(' Error: The center index should >0 and <=',i7)") ncenter
		end if
	
	else if (isel==22) then
		write(*,*) "Input index of the GTF, e.g. 13"
		read(*,*) i
		write(*,*) "Input the type you want to set, should be one of the following ones"
		write(*,*) "S,X,Y,Z,XX,YY,ZZ,XY,XZ,YZ,XXX,YYY,ZZZ,XXY,XXZ,YYZ,XYY,XZZ,YZZ,XYZ"
		write(*,*) "ZZZZ,YZZZ,YYZZ,YYYZ,YYYY,XZZZ,XYZZ,XYYZ,XYYY,XXZZ,XXYZ,XXYY,XXXZ,XXXY,XXXX"
		write(*,"(a)") " ZZZZZ,YZZZZ,YYZZZ,YYYZZ,YYYYZ,YYYYY,XZZZZ,XYZZZ,XYYZZ,XYYYZ,XYYYY,XXZZZ,XXYZZ,XXYYZ,XXYYY,XXXZZ,XXXYZ,XXXYY,XXXXZ,XXXXY,XXXXX"
		read(*,*) seltmpc
		do j=1,size(GTFtype2name)
			if (seltmpc==GTFtype2name(j)) then
				b(i)%type=j
				write(*,*) "Done!"
				exit
			end if
			if (j==20) write(*,*) "Error: Cannot recognize this type"
		end do
	
	else if (isel==23) then
		write(*,*) "Input index of the GTF, e.g. 13"
		read(*,*) i
		write(*,*) "Input exponent, e.g. 0.035"
		read(*,*) rexp
		b(i)%exp=rexp
        write(*,*) "Done!"
	
	else if (isel==24) then
		write(*,*) "Input index of the GTF, e.g. 13"
		read(*,*) iprm
		write(*,*) "Input orbital index, e.g. 12"
		read(*,*) imonum
		if (iprm<=nprims.and.iprm>0.and.imonum<=nmo.and.imonum>0) then
			write(*,*) "Input the coefficient, e.g. 2.83E-3"
			read(*,*) rcoeff
			CO(imonum,iprm)=rcoeff
            write(*,*) "Done!"
		else
			write(*,"(' Error: The index of the GTF or orbital exceeded valid range!')")
		end if
	
	else if (isel==25) then
		isetmode=1
		if (allocated(CObasa)) then
			write(*,*) "1 Set coefficients of some GTFs in some orbitals"
			write(*,*) "2 Set coefficients of some basis functions in some orbitals"
			read(*,*) isetmode
		end if
		if (isetmode==1) then
			allocate(idxsel(nprims)) 
            idxsel=1 !If idxsel(i)=1, then GTF i is selected
			write(*,"(a)") " You will be asked to input several conditions, the GTFs simultaneously satisfying them will be selected"
			write(*,*)
			write(*,*) "Input indices that the GTFs must be, e.g. 2,3,7-10"
            write(*,*) "Press ENTER button directly means ignoring this condition"
			read(*,"(a)") c2000tmp
            if (c2000tmp/=" ") then
				call str2arr(c2000tmp,ntmp)
				allocate(tmparrint(ntmp))
				call str2arr(c2000tmp,ntmp,tmparrint)
				do iGTF=1,nprims
					if (all(tmparrint/=iGTF)) idxsel(iGTF)=0
				end do
				deallocate(tmparrint)
            end if
			write(*,*) "Input indices of the atoms that the GTFs must be centered at, e.g. 3,9-13,18"
            write(*,*) "Press ENTER button directly means ignoring this condition"
			read(*,"(a)") c2000tmp
            if (c2000tmp/=" ") then
				call str2arr(c2000tmp,ntmp)
				allocate(tmparrint(ntmp))
				call str2arr(c2000tmp,ntmp,tmparrint)
				do iGTF=1,nprims
					if (all(tmparrint/=b(iGTF)%center)) idxsel(iGTF)=0
				end do
				deallocate(tmparrint)
            end if
			write(*,*) "Input the type of the GTFs must be, should be one of S,X,Y,Z,XX,XY"
			write(*,*) "You can also input one of S,P,D,F,G,H to select according to angular moment"
            write(*,"(a)") " Adding a negative sign means selecting all other types, for example, ""-P"" means selecting all GTFs other than P type"
            write(*,*) "Press ENTER button directly means ignoring this condition"
			read(*,"(a)") seltmpc
            if (seltmpc/=" ") then
				if (seltmpc(1:1)=='-') then
					seltmpc=seltmpc(2:)
					do iGTF=1,nprims
						if ( GTFtype2name(b(iGTF)%type)==trim(seltmpc).or.type2ang(b(iGTF)%type)==trim(seltmpc) ) idxsel(iGTF)=0
					end do
                else
					do iGTF=1,nprims
						if ( GTFtype2name(b(iGTF)%type)/=trim(seltmpc).and.type2ang(b(iGTF)%type)/=trim(seltmpc) ) idxsel(iGTF)=0
					end do
                end if
            end if
			write(*,"(' Coefficient of',i8,' GTFs have been selected')") count(idxsel==1)
			write(*,"(/,a)") " Input indices of the orbitals for which the coefficients of selected GTFs will be set, e.g. 2,9-13,18-25"
            write(*,*) "Press ENTER button directly selecting all orbitals"
			read(*,"(a)") c2000tmp
            if (c2000tmp/=" ") then
				call str2arr(c2000tmp,ntmp)
				allocate(tmparrint(ntmp))
				call str2arr(c2000tmp,ntmp,tmparrint)
            else
				ntmp=nmo
				allocate(tmparrint(ntmp))
                forall (iorb=1:nmo) tmparrint(iorb)=iorb
            end if
			write(*,*) "Input the expansion coefficient you want to set, e.g. 0.5"
			read(*,*) coval
            do iGTF=1,nprims
				if (idxsel(iGTF)==0) cycle
                CO(tmparrint(1:ntmp),iGTF)=coval
            end do
            deallocate(tmparrint,idxsel)
            write(*,*) "Done!"
		else if (isetmode==2) then
			allocate(idxsel(nbasis)) 
            idxsel=1 !If idxsel(i)=1, then basis function i is selected
			write(*,"(a)") " You will be asked to input several conditions, the basis function simultaneously satisfying them will be selected"
			write(*,*)
			write(*,*) "Input indices that the basis functions must be, e.g. 2,3,7-10"
            write(*,*) "Press ENTER button directly means ignoring this condition"
			read(*,"(a)") c2000tmp
            if (c2000tmp/=" ") then
				call str2arr(c2000tmp,ntmp)
				allocate(tmparrint(ntmp))
				call str2arr(c2000tmp,ntmp,tmparrint)
				do ibas=1,nbasis
					if (all(tmparrint/=ibas)) idxsel(ibas)=0
				end do
				deallocate(tmparrint)
            end if
			write(*,*) "Input indices of the atoms that the basis functions must be centered at"
            write(*,*) "For example, 3,9-13,18"
            write(*,*) "Press ENTER button directly means ignoring this condition"
			read(*,"(a)") c2000tmp
            if (c2000tmp/=" ") then
				call str2arr(c2000tmp,ntmp)
				allocate(tmparrint(ntmp))
				call str2arr(c2000tmp,ntmp,tmparrint)
				do ibas=1,nbasis
					if (all(tmparrint/=bascen(ibas))) idxsel(ibas)=0
				end do
				deallocate(tmparrint)
            end if
			write(*,*) "Input the type of the basis functions must be, should be one of S,X,Y,Z,XX,XY"
			write(*,*) "You can also input one of S,P,D,F,G,H to select according to angular moment"
            write(*,"(a)") " Adding a negative sign means selecting all other types, for example, ""-P"" means selecting all basis functions other than P type"
            write(*,*) "Press ENTER button directly means ignoring this condition"
			read(*,"(a)") seltmpc
            if (seltmpc/=" ") then
				if (seltmpc(1:1)=='-') then
					seltmpc=seltmpc(2:)
					do ibas=1,nbasis
						if ( GTFtype2name(bastype(ibas))==trim(seltmpc).or.type2ang(bastype(ibas))==trim(seltmpc) ) idxsel(ibas)=0
					end do
                else
					do ibas=1,nbasis
						if ( GTFtype2name(bastype(ibas))/=trim(seltmpc).and.type2ang(bastype(ibas))/=trim(seltmpc) ) idxsel(ibas)=0
					end do
                end if
            end if
			write(*,"(' Coefficient of',i8,' basis functions have been selected')") count(idxsel==1)
			write(*,"(/,a)") " Input indices of the orbitals for which the coefficients of selected basis functions will be set, e.g. 2,9-13,18-25"
            write(*,*) "Press ENTER button directly selecting all orbitals"
			read(*,"(a)") c2000tmp
            if (c2000tmp/=" ") then
				call str2arr(c2000tmp,ntmp)
				allocate(tmparrint(ntmp))
				call str2arr(c2000tmp,ntmp,tmparrint)
            else
				ntmp=nmo
				allocate(tmparrint(ntmp))
                forall (iorb=1:nmo) tmparrint(iorb)=iorb
            end if
			write(*,*) "Input the expansion coefficient you want to set, e.g. 0.5"
			read(*,*) coval
            do ibas=1,nbasis
				if (idxsel(ibas)==0) cycle
                do idxorb=1,ntmp
					iorb=tmparrint(idxorb)
                    if (iorb<=nbasis) then
						CObasa(ibas,iorb)=coval
                    else
						CObasb(ibas,iorb-nbasis)=coval
                    end if
                end do
            end do
            deallocate(tmparrint,idxsel)
            write(*,*) "Done!"
			imodwfn=1
		end if

	else if (isel==26) then
		do while(.true.)
			write(*,*)
			write(*,*) "Select the orbitals for which the occupation numbers are needed to be changed"
			write(*,*) "e.g. 2,4,13-16,20 means selecting orbitals 2,4,13,14,15,16,20"
			write(*,*) "Input 0 can select all orbitals, input q or 00 can return"
			read(*,"(a)") c1000tmp
			if (c1000tmp(1:1)=='q'.or.c1000tmp(1:2)=='00') exit
			if (c1000tmp(1:1)=='0') then
				numorbsel=nmo
				do i=1,nmo
					orbarr(i)=i
				end do
			else
				call str2arr(c1000tmp,numorbsel,orbarr)
				if ( any(orbarr(1:numorbsel)<1).or.any(orbarr(1:numorbsel)>nmo) ) then
					write(*,*) "Error: One or more orbital indices exceeded valid range!"
					cycle
				end if
			end if
            write(*,*)
			write(*,*) "Set occupation number to which value? e.g. 1.2"
			write(*,*) "Note:"
			write(*,"(a)") " You can also input for example ""+1.1"" ""-1.1"" ""*1.1"" ""/1.1"" to add, minus, multiply and divide the occupation numbers by 1.1"
			write(*,"(a)") " To recover the initial occupation numbers, input ""i"""
			write(*,"(a)") " To generate occupation state for calculating odd electron density, input ""odd"""
			read(*,"(a)") c1000tmp
			if (index(c1000tmp,"odd")/=0) then
                oddnum=0
				do iorb=1,numorbsel
					MOocc(orbarr(iorb))=min(2-MOocc(orbarr(iorb)),MOocc(orbarr(iorb)))
                    oddnum=oddnum+MOocc(orbarr(iorb))
				end do
				write(*,*) "Done!"
                write(*,"(a,f12.6)") " Sum of occupation numbers of selected orbitals:",oddnum
			else if (index(c1000tmp,"i")/=0) then
				MOocc(orbarr(1:numorbsel))=MOocc_org(orbarr(1:numorbsel))
				write(*,*) "The occupation numbers have been recovered"
			else if (c1000tmp(1:1)=='+'.or.c1000tmp(1:1)=='-'.or.c1000tmp(1:1)=='*'.or.c1000tmp(1:1)=='/') then
				read(c1000tmp(2:),*) tmpval
				if (c1000tmp(1:1)=='+') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))+tmpval
				if (c1000tmp(1:1)=='-') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))-tmpval
				if (c1000tmp(1:1)=='*') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))*tmpval
				if (c1000tmp(1:1)=='/') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))/tmpval
				write(*,*) "Done!"
			else
				read(c1000tmp,*) tmpval
				MOocc(orbarr(1:numorbsel))=tmpval
				write(*,*) "Done!"
			end if
			call updatenelec !Update the number of electrons
			imodwfn=1
			if (any(MOocc/=int(MOocc))) then
				if (wfntype==0) then
					wfntype=3 !RHF-> Restricted multiconfiguration wavefunction
					write(*,"(a)") " Note: Now the wavefunction is recognized as a restricted multiconfiguration wavefunction"
				else if (wfntype==1.or.wfntype==2) then !UHF/ROHF-> Unrestricted multiconfiguration wavefunction
					wfntype=4
					write(*,"(a)") " Note: Now the wavefunction is recognized as an unrestricted multiconfiguration wavefunction"
				end if
			end if
		end do
	
	else if (isel==27) then
		do while(.true.)
			write(*,*) "Set type for which range of orbitals?"
			write(*,*) "e.g. 2,4,13-16,20 means selecting orbitals 2,4,13,14,15,16,20"
			write(*,*) "Input 0 can select all orbitals, input q or 00 can return"
			read(*,"(a)") c1000tmp
			if (c1000tmp(1:1)=='q'.or.c1000tmp(1:2)=='00') exit
			if (c1000tmp(1:1)=='0') then
				numorbsel=nmo
				do i=1,nmo
					orbarr(i)=i
				end do
			else
				call str2arr(c1000tmp,numorbsel,orbarr)
				if ( any(orbarr(1:numorbsel)<1).or.any(orbarr(1:numorbsel)>nmo) ) then
					write(*,*) "Error: One or more orbital indices exceeded valid range!"
					cycle
				end if
			end if
			write(*,*) "Set to which type?  0=Alpha+Beta  1=Alpha  2=Beta"
			read(*,*) isettype
			MOtype(orbarr(1:numorbsel))=isettype
			!Recount alpha and beta electrons
			call updatenelec
			write(*,*) "Done!"
			!Update wavefunction type
			if (all(MOtype==0)) then
				if (all(MOocc==nint(MOocc))) then !All A+B orbital & integer occupation
					wfntype=0
					write(*,"(' Note: Now the wavefunction is recognized as a restricted closed-shell single-determinant wavefunction')")
				else !All A+B orbital & partial occupation
					wfntype=3
					write(*,"(' Note: Now the wavefunction is recognized as a restricted multiconfiguration wavefunction')")
				end if
			else
				if (any(MOocc/=nint(MOocc)).and.all(MOtype/=0)) then !Either A or B, and partial occupation
					wfntype=4
					write(*,"(' Note: Now the wavefunction is recognized as an unrestricted multiconfiguration wavefunction')")
				else if (all(MOocc==nint(MOocc)).and.all(MOtype/=0)) then !Integer occupation and either A or B
					wfntype=1
					write(*,"(' Note: Now the wavefunction is recognized as an unrestricted single-determinant wavefunction')")
				else if (all(MOocc==nint(MOocc)).and.any(MOtype==0).and.all(MOtype/=2)) then !Integer occupation and at least one orbital is A, and B is unexisted
					wfntype=2
					write(*,"(' Note: Now the wavefunction is recognized as a restricted open-shell wavefunction')")
				else
					write(*,"(' Warning: The type of present wavefunction cannot be identified! You need to reset orbital types')")
					write(*,*) "Press ENTER button to continue"
					read(*,*)
				end if
			end if
			imodwfn=1
		end do
		
	else if (isel==28) then
		do while(.true.)
			write(*,*) "Select the orbitals for which the energy are needed to be changed"
			write(*,*) "e.g. 2,4,13-16,20 means selecting orbitals 2,4,13,14,15,16,20"
			write(*,*) "Input 0 can select all orbitals, input q or 00 can return"
			read(*,"(a)") c1000tmp
			if (c1000tmp(1:1)=='q'.or.c1000tmp(1:2)=='00') exit
			if (c1000tmp(1:1)=='0') then
				numorbsel=nmo
				do i=1,nmo
					orbarr(i)=i
				end do
			else
				call str2arr(c1000tmp,numorbsel,orbarr)
				if ( any(orbarr(1:numorbsel)<1).or.any(orbarr(1:numorbsel)>nmo) ) then
					write(*,*) "Error: One or more orbital indices exceeded valid range!"
					cycle
				end if
			end if
			write(*,*) "0 Recover their initial orbital energies"
			write(*,*) "1 Set the orbital energies to a specific value"
			write(*,*) "2 Add a value to the orbital energies"
			write(*,*) "3 Minus a value from the orbital energies"
			write(*,*) "4 Multiply the orbital energies by a value"
			write(*,*) "5 Divide the orbital energies by a value"
			read(*,*) iselop
			if (iselop==0) then
				MOene(orbarr(1:numorbsel))=MOene_org(orbarr(1:numorbsel))
			else if (iselop==1) then
				write(*,*) "Input the energy in eV, e.g. -3.6"
				read(*,*) tmpval
				MOene(orbarr(1:numorbsel))=tmpval/au2eV
			else if (iselop==2) then
				write(*,*) "Input the value in eV, e.g. 0.6"
				read(*,*) tmpval
				MOene(orbarr(1:numorbsel))=MOene(orbarr(1:numorbsel))+tmpval/au2eV
			else if (iselop==3) then
				write(*,*) "Input the value in eV, e.g. 0.6"
				read(*,*) tmpval
				MOene(orbarr(1:numorbsel))=MOene(orbarr(1:numorbsel))-tmpval/au2eV
			else if (iselop==4) then
				write(*,*) "Input the factor, e.g. 0.9"
				read(*,*) tmpval
				MOene(orbarr(1:numorbsel))=MOene(orbarr(1:numorbsel))*tmpval
			else if (iselop==5) then
				write(*,*) "Input the factor, e.g. 1.2"
				read(*,*) tmpval
				MOene(orbarr(1:numorbsel))=MOene(orbarr(1:numorbsel))/tmpval
			end if
			write(*,*) "Done!"
		end do
	
	else if (isel==29) then
		write(*,"(a)") " Input indices of two orbitals, e.g. 3,8, then all of their information (energy, occupation, coefficients) will be exchanged"
        read(*,*) iorb,jorb
        if (MOtype(iorb)==MOtype(jorb)) then
			enetmp=MOene(iorb)
			MOene(iorb)=MOene(jorb)
			MOene(jorb)=enetmp
			occtmp=MOocc(iorb)
			MOocc(iorb)=MOocc(jorb)
			MOocc(jorb)=enetmp
			allocate(tmparr(nprims))
			tmparr(:)=CO(iorb,:)
			CO(iorb,:)=CO(jorb,:)
			CO(jorb,:)=tmparr(:)
			deallocate(tmparr)
			if (allocated(CObasa)) then
				allocate(tmparr(nbasis))
				if (iorb<=nbasis) then
					tmparr(:)=CObasa(:,iorb)
					CObasa(:,iorb)=CObasa(:,jorb)
					CObasa(:,jorb)=tmparr(:)
                else
					tmparr(:)=CObasb(:,iorb-nbasis)
					CObasb(:,iorb-nbasis)=CObasb(:,jorb-nbasis)
					CObasb(:,jorb-nbasis)=tmparr(:)
                end if
				deallocate(tmparr)
			end if
			imodwfn=1
			write(*,*) "Done!"
        else
			write(*,*) "Error: In order to alter information, the two orbitals must have the same type!"
            write(*,*) "Press ENTER button to continue"
            read(*,*)
        end if
	
	else if (isel==30) then
		write(*,*) "1 Exchange orbital energies (in eV) with occupation numbers"
		write(*,*) "2 Exchange orbital energies (in Hartree) with occupation numbers"
        write(*,*) "Note: For EDDB user, you should choose 1"
        read(*,*) isel2
		allocate(tmparr(nmo))
        tmparr=MOocc
        if (isel2==1) then
			MOocc=MOene*au2eV
        else if (isel2==2) then
			MOocc=MOene
        end if
        MOene=tmparr
        deallocate(tmparr)
        write(*,*) "Orbital energies and occupation numbers have been exchanged with each other"
		if (any(MOocc/=int(MOocc))) then
			if (wfntype==0) then
				wfntype=3 !RHF-> Restricted multiconfiguration wavefunction
				write(*,"(a)") " Note: Now the wavefunction is recognized as a restricted multiconfiguration wavefunction"
			else if (wfntype==1.or.wfntype==2) then !UHF/ROHF-> Unrestricted multiconfiguration wavefunction
				wfntype=4
				write(*,"(a)") " Note: Now the wavefunction is recognized as an unrestricted multiconfiguration wavefunction"
			end if
		end if
        imodwfn=1
        call updatenelec
    
	else if (isel==31) then
		write(*,*) "Input X,Y,Z of translation vector (e.g. 3.2,1.0,0)"
		read(*,*) pbctransx,pbctransy,pbctransz
		write(*,*) "You inputted coordinates are in which unit?  1: Bohr  2: Angstrom"
		read(*,*) iunit
		if (iunit==2) then
			pbctransx=pbctransx/b2a
			pbctransy=pbctransy/b2a
			pbctransz=pbctransz/b2a
		end if
		do i=1,ncenter
			a(i)%x=a(i)%x+pbctransx
			a(i)%y=a(i)%y+pbctransy
			a(i)%z=a(i)%z+pbctransz
		end do
		imodwfn=1
	
	else if (isel==32) then
		write(*,*) "Input X,Y,Z of translation vector (e.g. 3.2,1.0,0)"
		read(*,*) pbctransx,pbctransy,pbctransz
		write(*,*) "You inputted coordinates are in which unit?  1: Bohr  2: Angstrom"
		read(*,*) iunit
		if (iunit==2) then
			pbctransx=pbctransx/b2a
			pbctransy=pbctransy/b2a
			pbctransz=pbctransz/b2a
		end if
		write(*,*) "Duplicate system how many times? e.g. 3"
		read(*,*) numdup
		!_tmp is for backing up current information
        if (allocated(a_tmp)) deallocate(a_tmp)
		allocate(a_tmp(ncenter))
		allocate(b_tmp(nprims))
		allocate(CO_tmp(nmo,nprims))
		a_tmp=a
		b_tmp=b
		CO_tmp=CO
		deallocate(a,b,CO)
		nprims_tmp=nprims
		ncenter_tmp=ncenter
		nprims=nprims*(numdup+1)
		ncenter=ncenter*(numdup+1)
		nelec=nelec*(numdup+1)
		naelec=naelec*(numdup+1)
		nbelec=nbelec*(numdup+1)
		allocate(a(ncenter))
		allocate(b(nprims))
		allocate(CO(nmo,nprims))
		do idup=0,numdup
			a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))=a_tmp(1:ncenter_tmp)
			a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))%x=a_tmp(1:ncenter_tmp)%x+pbctransx*idup
			a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))%y=a_tmp(1:ncenter_tmp)%y+pbctransy*idup
			a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))%z=a_tmp(1:ncenter_tmp)%z+pbctransz*idup
			b(nprims_tmp*idup+1:nprims_tmp*(idup+1))=b_tmp(1:nprims_tmp)
			b(nprims_tmp*idup+1:nprims_tmp*(idup+1))%center=b_tmp(1:nprims_tmp)%center+ncenter_tmp*idup
			CO(:,nprims_tmp*idup+1:nprims_tmp*(idup+1))=CO_tmp(:,1:nprims_tmp) !Notice that the orbitals do not satisify normalization condition any more, and the orbital occupation number will be artifical
		end do
		deallocate(a_tmp,b_tmp,CO_tmp)
		imodwfn=1
		call gendistmat !The number of atoms have changed, so we must update distance matrix
	
	else if (isel==33) then
		write(*,*) "Rotate which orbital? (Input 0 to rotate all orbitals)"
		read(*,*) iorb
		if (iorb/=0) then
			call orbcoeffrotate(iorb)
		else if (iorb==0) then
            write(*,*) "Please wait..."
			do imo=1,nmo
				call orbcoeffrotate(imo)
			end do
			write(*,*) "Also rotate atomic coordinates? (y/n)"
			read(*,*) selectyn
			if (selectyn=='y'.or.selectyn=='Y') then
				do iatm=1,ncenter
					tmpval=a(iatm)%x
					a(iatm)%x=a(iatm)%z
					a(iatm)%z=a(iatm)%y
					a(iatm)%y=tmpval
				end do
			end if
		end if
		write(*,*) "Done!"
	
	else if (isel==34) then
		call getninnerele(ninnerele,1)
		nelec=nelec-ninnerele
		naelec=naelec-ninnerele/2
		nbelec=nbelec-ninnerele/2
		if (wfntype==1.or.wfntype==4) then !UHF and U-post-HF wfn
			MOocc(1:ninnerele/2)=0D0
			do j=1,nmo !Where the first beta orbital appear now
				if (motype(j)==2) exit
			end do
			MOocc(j:j+ninnerele/2-1)=0D0
			write(*,"(' The occupation of',i7,' lowest energy orbitals have been set to zero')") ninnerele
		else if (wfntype==0.or.wfntype==2.or.wfntype==3) then !restricted(=0) or RO(=2) or post-R(=3) wavefunction
			MOocc(1:ninnerele/2)=0D0
			write(*,"(' The occupation of',i7,' lowest energy orbitals have been set to zero')") ninnerele/2
		end if
		if (wfntype==3.or.wfntype==4) write(*,"(' Warning: Discarding inner orbitals for multiconfiguration wavefunction will lead to unexpected result!')") 
		imodwfn=1
	
	else if (isel==35) then
		call selMO_IRREP
	
	else if (isel==36) then
		write(*,*) "Input index of the orbitals, e.g. 2,3,7-10"
		read(*,"(a)") c1000tmp
		call str2arr(c1000tmp,ntmp,orbarr)
		do idxtmp=1,ntmp
			idx=orbarr(idxtmp)
			CO(idx,:)=-CO(idx,:)
			if (allocated(CObasa)) then
				if (idx<=nbasis) then
					CObasa(:,idx)=-CObasa(:,idx)
				else
					CObasb(:,idx-nbasis)=-CObasb(:,idx-nbasis)
				end if
			end if
		end do
		write(*,*) "Done!"
		imodwfn=1
	
	else if (isel==37) then
		if (wfntype==1.or.wfntype==4) then
			write(*,*) "Error: The current wavefunction is already unrestricted!"
            write(*,*) "Press ENTER button to return"
            read(*,*)
            cycle
        end if
        !Back up current wavefunction
		if (allocated(CO_back)) deallocate(CO_back,MOene_back,MOocc_back)
		allocate(CO_back(nmo,nprims),MOene_back(nmo),MOocc_back(nmo))
		nmo_back=nmo
		CO_back=CO
		MOene_back=MOene
		MOocc_back=MOocc
        !Reallocate arrays with new size
        deallocate(CO,MOene,MOocc,MOtype)
        nmo=nmo_back*2
        allocate(CO(nmo,nprims),MOene(nmo),MOocc(nmo),MOtype(nmo))
        if (wfntype==0.or.wfntype==3) then !RHF, R-post-HF
			do imo=1,nmo_back
				MOene(imo)=MOene_back(imo)
				MOene(imo+nmo_back)=MOene_back(imo)
				MOocc(imo)=MOocc_back(imo)/2
				MOocc(imo+nmo_back)=MOocc_back(imo)/2
                MOtype(imo)=1
                MOtype(imo+nmo_back)=2
                CO(imo,:)=CO_back(imo,:)
                CO(imo+nmo_back,:)=CO_back(imo,:)
            end do
            if (wfntype==0) wfntype=1
            if (wfntype==3) wfntype=4
        else if (wfntype==2) then !ROHF
			MOocc=0
			do imo=1,nmo_back
				MOene(imo)=MOene_back(imo)
				MOene(imo+nmo_back)=MOene_back(imo)
                if (MOocc_back(imo)==2) then
					MOocc(imo)=1
					MOocc(imo+nmo_back)=1
                else if (MOocc_back(imo)==1) then
					MOocc(imo)=1
                end if
                MOtype(imo)=1
                MOtype(imo+nmo_back)=2
                CO(imo,:)=CO_back(imo,:)
                CO(imo+nmo_back,:)=CO_back(imo,:)
            end do
            wfntype=1
        end if
        deallocate(CO_back,MOene_back,MOocc_back)
        if (allocated(CObasa)) then
			allocate(CObasb(nbasis,nbasis))
            CObasb=CObasa
        end if
        imodwfn=1
        write(*,*) "Finished!"
	end if
    
end do
end subroutine




!!---------------- Select MOs according to irreducible representation
subroutine selMO_IRREP
use defvar
use util
character symlab(nmo)*4,c2000tmp*2000,symstat(nmo)*9 !Allocate the array lengths as upper limit
integer tmparr(nmo),symNorb(nmo) !Allocate the array lengths as upper limit
if (wfntype/=0.and.wfntype/=1) then
	write(*,"(a)") " Error: This function only works for RHF or UHF wavefunctions (or the DFT counterparts)"
	return
end if
nsym=0
do imo=1,nmo
	if (MOocc_org(imo)==0D0) cycle
	if (all(symlab(1:nsym)/=MOsym(imo))) then
		nsym=nsym+1
		symlab(nsym)=MOsym(imo)
	end if
end do
symNorb=0
do imo=1,nmo
	if (MOocc_org(imo)==0D0) cycle
	do isym=1,nsym
		if (MOsym(imo)==symlab(isym)) symNorb(isym)=symNorb(isym)+1
	end do
end do
symstat="Normal"
MOocc=MOocc_org
if (imodwfn==1) write(*,*) "Note: Original occupation status has been recovered"
write(*,*) "Note: Only the orbitals that originally occupied are taken into account here"
do while(.true.)
	write(*,*)
	write(*,*) "Information of various irreducible representations:"
	do isym=1,nsym
		write(*,"(i5,'  Sym: ',a,'  N_orb:',i5,'    Status: ',a)") isym,symlab(isym),symNorb(isym),symstat(isym)
	end do
	write(*,*)
	write(*,*) "0 Save and return"
	write(*,*) "1 Discard specific irreducible representations"
	write(*,*) "2 Recover original status"
	write(*,*) "3 Reverse status"
	read(*,*) isel
	
	if (isel==0) then
		call updatenelec
		imodwfn=1
		write(*,*) "The current orbital occupation status has been saved"
		write(*,*) "Updating density matrix..."
		call genP
		write(*,*) "Density matrix has been updated"
		exit
	else if (isel==2) then
		MOocc=MOocc_org
		symstat="Normal"
	else if (isel==1.or.isel==3) then
		if (isel==1) then
			write(*,*) "Input the index of the irreducible representations to be discarded, e.g. 1,3-5"
			read(*,"(a)") c2000tmp
			call str2arr(c2000tmp,nsymsel,tmparr)
			do isym=1,nsymsel
				symstat(tmparr(isym))="Discarded"
			end do
		else if (isel==3) then
			do isym=1,nsym
				if (symstat(isym)=="Normal") then
					symstat(isym)="Discarded"
				else
					symstat(isym)="Normal"
				end if
			end do
		end if
		do imo=1,nmo
			if (MOocc_org(imo)==0D0) cycle
			do isym=1,nsym
				if (MOsym(imo)==symlab(isym)) then
					if (symstat(isym)=="Normal") MOocc(imo)=MOocc_org(imo)
					if (symstat(isym)=="Discarded") MOocc(imo)=0D0
					exit
				end if
			end do
		end do
	end if
	write(*,*) "Done!"
end do
end subroutine



!!---------- Update the number of electrons
subroutine updatenelec
use defvar
integer imo
nelec=0
naelec=0
nbelec=0
do imo=1,nmo
	if (MOtype(imo)==0) then
		naelec=naelec+MOocc(imo)/2D0
		nbelec=nbelec+MOocc(imo)/2D0
	else if (MOtype(imo)==1) then
		naelec=naelec+MOocc(imo)
	else if (MOtype(imo)==2) then
		nbelec=nbelec+MOocc(imo)
	end if
end do
nelec=naelec+nbelec
end subroutine



!!-------- Generate nelec, naelec, nbelec by guessing when wavefunction information is not available
!Commonly invoked when reading file only contains geometry information
!The system is always regarded as neutral. If number of electrons is even, then naelec=nbelec=nelec/2. Else odd, naelec=nbelec+1
subroutine guessnelec
use defvar
nelec=sum(a%charge)
if (mod(nint(nelec),2)==0) then
    naelec=nint(nelec)/2
    nbelec=naelec
else
    nbelec=(nint(nelec)-1)/2
    naelec=nbelec+1
end if

end subroutine
			
			

!!-------- Check if present wavefunction is sanity, i.e. all orbital satisfies normalization condition
subroutine wfnsanity
use defvar
implicit real*8 (a-h,o-z)
real*8 GTFSmat(nprims*(nprims+1)/2)

rmaxdev=0
rmaxdevint=0
if (ifPBC==0) then
    call genGTFSmat(GTFSmat,nprims*(nprims+1)/2)
    do imo=1,nmo
	    tmp=0
	    !$OMP parallel shared(tmp) private(iGTF,jGTF,tmpprivate) num_threads(nthreads)
	    tmpprivate=0
	    !$OMP do schedule(dynamic)
	    do iGTF=1,nprims
		    do jGTF=iGTF+1,nprims
			    tmpprivate=tmpprivate+2*CO(imo,iGTF)*CO(imo,jGTF)*GTFSmat(jGTF*(jGTF-1)/2+iGTF)
		    end do
		    tmpprivate=tmpprivate+CO(imo,iGTF)**2*GTFSmat(iGTF*(iGTF-1)/2+iGTF)
	    end do
	    !$OMP END DO
	    !$OMP CRITICAL
	    tmp=tmp+tmpprivate
	    !$OMP END CRITICAL
	    !$OMP END PARALLEL
	    write(*,"(' Orbital',i7,', Occ:',f8.4,'   Value:',f16.10)") imo,MOocc(imo),tmp
	    tmpt=abs(tmp-1)
	    if (tmpt>rmaxdev) rmaxdev=tmpt
	    tmpt=abs(tmp-nint(tmp))
	    if (tmpt>rmaxdevint) rmaxdevint=tmpt
    end do

else !Examine using basis function information
    call ask_Sbas_PBC
    do imo=1,nbasis
	    tmp=0
        if (all(CObasa(:,imo)==0)) cycle
        do ibas=1,nbasis
            do jbas=1,nbasis
                tmp=tmp+CObasa(ibas,imo)*CObasa(jbas,imo)*Sbas(ibas,jbas)
            end do
        end do
	    write(*,"(' Orbital',i7,', Occ:',f8.4,'   Value:',f16.10)") imo,MOocc(imo),tmp
	    tmpt=abs(tmp-1)
	    if (tmpt>rmaxdev) rmaxdev=tmpt
	    tmpt=abs(tmp-nint(tmp))
	    if (tmpt>rmaxdevint) rmaxdevint=tmpt
    end do
end if
write(*,"(' Maximum deviation to 1:',f16.10)") rmaxdev
write(*,"(' Maximum deviation to integer:',f16.10)") rmaxdevint

write(*,*) "Press ENTER button to continue"
read(*,*)
end subroutine




!!---------- Return normalization coefficient for specific type of Cartesian type GTF, see Levine 5ed p487
!The meaning of itype is defined in GTFtype2name
real*8 function normgau(itype,exp)
use defvar
use util
implicit real*8 (a-h,o-z)
real*8 exp
ix=type2ix(itype)
iy=type2iy(itype)
iz=type2iz(itype)
normgau=(2*exp/pi)**0.75D0*dsqrt( (8*exp)**(ix+iy+iz)*ft(ix)*ft(iy)*ft(iz)/(ft(2*ix)*ft(2*iy)*ft(2*iz)) )
end function


!---- Renormalizing shells (modifying contraction coefficients in basis shells) for Molden input file
!Various basis functions in a shell may have different normalization factor, however, since we are dealing with shell, 
!we can use any specific basis function to derive normalization factor and fix contraction coefficients
!Here we only consider (Lval,0,0) or (0,Lval,0) or (0,0,Lval) type, which have identical normalization factor
!Lval=0/1/2/3/4/5 corresponds to s/p/d/f/g/h, exp and con are exponents and contraction coefficients in the basis shell, respectively
subroutine renormmoldengau(nlen,Lval,exp,con)
implicit real*8 (a-h,o-z)
integer nlen,Lval
real*8 exp(nlen),con(nlen),ctmp(nlen)
pi=acos(-1D0)
call genn1n2nf(Lval,n1,n2,nf)
fc=2D0**n1/(pi**3*nf)
do i=1,nlen
	prmnormfac=sqrt(sqrt(fc*exp(i)**n2)) !Normalization factor of (Lval,0,0) type of GTF
	ctmp(i)=con(i)*prmnormfac
end do
facnorm=0D0 !Calculate <Lval,0,0|Lval,0,0> overlap integral for (Lval,0,0) type of contracted basis function
do i=1,nlen
	do j=1,i
	  expavg=(exp(i)+exp(j))/2D0
	  facadd=ctmp(i)*ctmp(j)/sqrt(fc*expavg**n2)
	  if (i/=j) facadd=facadd*2
	  facnorm=facnorm+facadd
	end do
end do
if (facnorm>1D-10) facnorm=1/sqrt(facnorm)
con=con*facnorm !Fix contraction coefficients
end subroutine
!!---- Produce normalization factor of (Lval,0,0) type of GTF with angular moment of Lval, &
!used to fix the contraction coefficient problem of the Molden input file generated by ORCA
!This routine can also be replaced by the slower routine "normgau", e.g. renormgau_ORCA(0.8D0,3) = normgau(21,0.8D0)
!cf. fnorm_lmn of m2a
real*8 function renormgau_ORCA(exp,Lval)
real*8 exp,pi
integer Lval,n1,n2,nf
pi=acos(-1D0)
call genn1n2nf(Lval,n1,n2,nf)
renormgau_ORCA=dsqrt(dsqrt(2**n1*exp**n2/(pi**3*nf))) !Get norm, the norm^4 = 2^n1 * a^n2 / (pi^3 * nf)
end function
!!--- Generate n1,n2,f for a given angular moment (up to h), the data only correspond to (Lval,0,0)=(0,Lval,0)=(0,0,Lval) case
!n1=3+4*(l+m+n)
!n2=3+2*(l+m+n)
!nf=[(2l-1)!!(2m-1)!!(2n-1)!!]^2  <--- e.g. f shell, putting XXXX, i.e. (4,0,0) into it returns 11025. Here 7!! means 7*5*3*1
subroutine genn1n2nf(Lval,n1,n2,nf)
integer Lval,n1,n2,nf
if (Lval==0) then
	n1=3
	n2=3
	nf=1
else if (Lval==1) then
	n1=7
	n2=5
	nf=1
else if (Lval==2) then
	n1=11
	n2=7
	nf=9
else if (Lval==3) then
	n1=15
	n2=9
	nf=225
else if (Lval==4) then
	n1=19
	n2=11
	nf=11025
else if (Lval==5) then
	n1=23
	n2=13
	nf=893025
end if
end subroutine



!!----- Use Lowdin orthogonalization method to transform density matrix and coefficient matrix &
!!----- to orthonormal basis, meanwhile update Sbas to identity matrix
!See Szabo p143 for details
subroutine symmortho
use defvar
use util
implicit real*8 (a-h,o-z)
real*8 Umat(nbasis,nbasis),svalvec(nbasis),Xmat(nbasis,nbasis) !workvec(3*nbasis-1)
! call DSYEV('V','U',nbasis,sbas,nbasis,svalvec,workvec,3*nbasis-1,ierror) !lapack, the resultant sbas is eigenvector matrix
call diagsymat(Sbas,Umat,svalvec,ierror)
if (ierror/=0) write(*,*) "Error: Diagonalization of overlap matrix failed!"
!Now Sbas is already diagonalized
forall (i=1:nbasis) Sbas(i,i)=dsqrt(svalvec(i)) !Use Sbas as temporary matrix here
Xmat=matmul(matmul(Umat,Sbas),transpose(Umat)) !Then Xmat is S^0.5
Ptot=matmul(matmul(Xmat,Ptot),Xmat)
if (allocated(Palpha)) then
	Palpha=matmul(matmul(Xmat,Palpha),Xmat)
	Pbeta=Ptot-Palpha
end if
Cobasa=matmul(Xmat,Cobasa)
if (allocated(Cobasb)) Cobasb=matmul(Xmat,Cobasb)
forall(i=1:nbasis) Sbas(i,i)=1D0 !Reconstruct overlap matrix in orthonormal basis function
end subroutine

!!--- Input overlap matrix and return Lowdin orthogonalization transformation matrix Xmat=S^0.5 and Xmatinv=S^-0.5
! Smatin is inputted overlap matrix, which will not be modified
subroutine symmorthomat(Smatin,Xmat,Xmatinv)
use defvar
use util
real*8 Umat(nbasis,nbasis),svalvec(nbasis),Smatin(nbasis,nbasis),Smat(nbasis,nbasis),Xmat(nbasis,nbasis),Xmatinv(nbasis,nbasis)
Smat=Smatin
call diagsymat(Smat,Umat,svalvec,ierror)
if (ierror/=0) write(*,*) "Error: Diagonalization of overlap matrix is fail!"
forall (i=1:nbasis) Smat(i,i)=dsqrt(svalvec(i))
Xmat=matmul(matmul(Umat,Smat),transpose(Umat))
forall (i=1:nbasis) Smat(i,i)=1D0/Smat(i,i)
Xmatinv=matmul(matmul(Umat,Smat),transpose(Umat))
end subroutine


!!!------------------ Generate distance matrix between atoms (in Bohr)
!To avoid complicate thing, PBC is not taken into account
subroutine gendistmat
use defvar
implicit real*8 (a-h,o-z)
if (allocated(distmat)) deallocate(distmat)
allocate(distmat(ncenter,ncenter))
distmat=0.0D0
!if (ifPBC==0) then
	!$OMP PARALLEL DO SHARED(distmat) PRIVATE(i,j,tmp) schedule(dynamic) NUM_THREADS(nthreads)
	do i=1,ncenter
		do j=i+1,ncenter
			tmp=dsqrt((a(i)%x-a(j)%x)**2+(a(i)%y-a(j)%y)**2+(a(i)%z-a(j)%z)**2)
			distmat(i,j)=tmp
			distmat(j,i)=tmp
		end do
	end do
	!$OMP END PARALLEL DO
!else
!	!$OMP PARALLEL DO SHARED(distmat) PRIVATE(i,j) schedule(dynamic) NUM_THREADS(nthreads)
!	do i=1,ncenter
!		do j=i+1,ncenter
!			call nearest_atmdistxyz(i,j,distmat(i,j),atmx,atmy,atmz)
!			distmat(j,i)=distmat(i,j)
!		end do
!	end do
!	!$OMP END PARALLEL DO
!end if

end subroutine


!!!------------ Return the number of inner-core orbitals
subroutine getninnerele(ninnerele,info)
use defvar
integer ninnerele,info
ninnerele=0
do i=1,ncenter
	if (int(a(i)%charge)/=a(i)%index) then
		if (info==1) write(*,"(' Note: Atom',i5,' is not taken into account since it utilizes pseudopotential')") i
		cycle
	end if
	if (a(i)%index>2.and.a(i)%index<=10) ninnerele=ninnerele+2
	if (a(i)%index>10.and.a(i)%index<=18) ninnerele=ninnerele+10
	if (a(i)%index>18.and.a(i)%index<=36) ninnerele=ninnerele+18
	if (a(i)%index>36.and.a(i)%index<=54) ninnerele=ninnerele+36
	if (a(i)%index>54.and.a(i)%index<=86) ninnerele=ninnerele+54
	if (a(i)%index>86) ninnerele=ninnerele+86
end do
end subroutine


!!!------------------------- Swap two GTF
subroutine swapGTF(i,j,swaptype)
use defvar
integer n,i,j
character swaptype*3
type(primtype) tempb !For exchanging basis functions' order
if (swaptype=="all") then
	tempb=b(i)
	b(i)=b(j)
	b(j)=tempb
else if (swaptype=="cen") then
	tempb%center=b(i)%center
	b(i)%center=b(j)%center
	b(j)%center=tempb%center
else if (swaptype=="typ") then
	tempb%type=b(i)%type
	b(i)%type=b(j)%type
	b(j)%type=tempb%type
else if (swaptype=="exp") then
	tempb%exp=b(i)%exp
	b(i)%exp=b(j)%exp
	b(j)%exp=tempb%exp
end if
if (swaptype=="all".or.swaptype=="MO ") then
	do n=1,nmo
		temp=CO(n,i)
		CO(n,i)=CO(n,j)
		CO(n,j)=temp
	end do
end if
end subroutine


!!!---- Rotate(exchange) GTF and basis function coefficients within all shell in different direction of specific orbital
! If using this three times, namely XYZ->ZXY->YZX->XYZ, the original coefficient will be recovered
! In detail, for example, for d-type will lead to such coefficient exchange: XX to YY, YY to ZZ, ZZ to XX, XY to YZ, XZ to XY, YZ to XZ
! exchange only involve the GTFs/basis func. in the same shell
subroutine orbcoeffrotate(orb) !orb=Rotate which orbital (from 1 to nmo)
use defvar
implicit real*8 (a-h,o-z)
integer orb
real*8 COorborg(nprims) !For backing up origin CO
real*8 CObasa_tmp(nbasis) !For backing up origin CObasa
real*8 CObasb_tmp(nbasis) !For backing up origin CObasb
COorborg(:)=CO(orb,:)
do i=1,nprims
	ixtmp=type2iz(b(i)%type)
	iytmp=type2ix(b(i)%type)
	iztmp=type2iy(b(i)%type)
	do j=1,nprims
		if (type2ix(b(j)%type)==ixtmp.and.type2iy(b(j)%type)==iytmp.and.&
		type2iz(b(j)%type)==iztmp.and.b(j)%exp==b(i)%exp.and.b(j)%center==b(i)%center) CO(orb,j)=COorborg(i)
	end do
end do
if (allocated(CObasa)) then
    if (orb<=nbasis) then
	    CObasa_tmp=CObasa(:,orb)
	else
        CObasb_tmp=CObasb(:,orb-nbasis)
    end if
	do iatm=1,ncenter
		do ibas=basstart(iatm),basend(iatm)
			ityp=bastype(ibas)
			ixtmp=type2iz(ityp)
			iytmp=type2ix(ityp)
			iztmp=type2iy(ityp)
			do jbas=basstart(iatm),basend(iatm)
				jtyp=bastype(jbas)
				if (type2ix(jtyp)==ixtmp.and.type2iy(jtyp)==iytmp.and.type2iz(jtyp)==iztmp.and.basshell(ibas)==basshell(jbas)) then
                    if (orb<=nbasis) then
					    CObasa(jbas,orb)=CObasa_tmp(ibas)
					else
                        CObasb(jbas,orb-nbasis)=CObasb_tmp(ibas)
                    end if
				end if
			end do
		end do
	end do
end if
end subroutine



!!!------ Output molecular formula
subroutine showformula
use defvar
implicit real*8 (a-h,o-z)
character tmp*6
write(*,"(' Formula: ')",advance="no")
do i=0,nelesupp
	n=0
	do iatm=1,ncenter
		if (a(iatm)%index==i) n=n+1
	end do
	write(tmp,"(i6)") n
	if (n/=0) write(*,"(a,a,' ')",advance="no") trim(ind2name(i)),trim(adjustl(tmp))
end do
write(*,*)
end subroutine



!!!----------- Resize number of orbitals of CO, MOene, MOtype, MOocc to "newnmo", also resize number of GTFs of CO to "newnprims"
subroutine resizebynmo(newnmo,newnprims)
use defvar
implicit real*8 (a-h,o-z)
real*8,allocatable :: CO_bk(:,:),MOene_bk(:),MOocc_bk(:)
integer,allocatable :: MOtype_bk(:)
integer newnmo,oldnmo,newnprims,oldnprims
oldnmo=size(CO,1)
oldnprims=size(CO,2)
allocate(CO_bk(oldnmo,oldnprims),MOene_bk(oldnmo),MOocc_bk(oldnmo),MOtype_bk(oldnmo))
CO_bk=CO
MOene_bk=MOene
MOocc_bk=MOocc
MOtype_bk=MOtype
deallocate(CO,MOene,MOocc,MOtype)
allocate(CO(newnmo,newnprims),MOene(newnmo),MOocc(newnmo),MOtype(newnmo))
if (newnmo>=oldnmo) then !Enlarge array size, don't forget to fill the gap afterwards
	if (newnprims>=oldnprims) CO(1:oldnmo,1:oldnprims)=CO_bk(:,:)
	if (newnprims<oldnprims) CO(1:oldnmo,:)=CO_bk(:,1:newnprims)
	MOene(1:oldnmo)=MOene_bk(:)
	MOocc(1:oldnmo)=MOocc_bk(:)
	MOtype(1:oldnmo)=MOtype_bk(:)
else if (newnmo<oldnmo) then !Reduce array size
	if (newnprims>=oldnprims) CO(:,1:oldnprims)=CO_bk(1:newnmo,:)
	if (newnprims<oldnprims) CO(:,:)=CO_bk(1:newnmo,1:newnprims)
	MOene(:)=MOene_bk(1:newnmo)
	MOocc(:)=MOocc_bk(1:newnmo)
	MOtype(:)=MOtype_bk(1:newnmo)
end if
deallocate(CO_bk,MOene_bk,MOocc_bk,MOtype_bk)
end subroutine


!!!------------ Generate gjf of atoms in molecule, and invoke Gaussian to get .wfn, then input them into custom list
subroutine setpromol
use defvar
use util
implicit real*8 (a-h,o-z)
integer :: itype=0
character(len=2) typename(100),nametmp
character basisset*80,tmpdir*80,c80tmp*80,outwfnname*80
logical alivegauout,alivewfntmp,aliveatomwfn
if (isys==1) call delfile("gxx.* fort.6 Gau*.inp") !Clean Gaussian scratch files in current folder

!The only difference between c80tmp and tmpdir is that the latter has \ or / separator at the end
if (iwfntmptype==1) then
	if (isys==1) tmpdir="wfntmp\"
	if (isys==2) tmpdir="wfntmp/"
	c80tmp="wfntmp"
	inquire(directory="wfntmp",exist=alivewfntmp)
	if (isys==1.and.alivewfntmp) then !Delete old wfntmp folder
		write(*,*) "Running: rmdir /S /Q wfntmp"
		call system("rmdir /S /Q wfntmp")
	else if (isys==2.and.alivewfntmp) then
		write(*,*) "Running: rm -rf wfntmp"
		call system("rm -rf wfntmp")
	end if
else if (iwfntmptype==2) then
	do i=1,9999 !Find a proper name of temporary folder
		write(c80tmp,"('wfntmp',i4.4)") i
		inquire(directory=c80tmp,exist=alivewfntmp)
		if (.not.alivewfntmp) exit
	end do
	if (isys==1) write(tmpdir,"('wfntmp',i4.4,'\')") i
	if (isys==2) write(tmpdir,"('wfntmp',i4.4,'/')") i
end if
write(*,*) "Running: mkdir "//trim(c80tmp) !Build new temporary folder
call system("mkdir "//trim(c80tmp))
inquire(directory="atomwfn",exist=aliveatomwfn)
if (isys==1.and.aliveatomwfn) then
	write(*,*) "Running: copy atomwfn\*.wfn "//trim(tmpdir)
	call system("copy atomwfn\*.wfn "//trim(tmpdir))
else if (isys==2.and.aliveatomwfn) then
	write(*,*) "Running: cp atomwfn/*.wfn "//trim(tmpdir)
	call system("cp atomwfn/*.wfn "//trim(tmpdir))
end if

noatmwfn=0 !Check if the atomic wfn file have pre-stored in atomwfn folder, if not, invoke gaussian to calc it
do i=1,nfragatm
	if (isys==1) inquire(file="atomwfn\"//a(fragatm(i))%name//".wfn",exist=alive)
	if (isys==2) inquire(file="atomwfn/"//a(fragatm(i))%name//".wfn",exist=alive)
	if (.not.alive) then
		noatmwfn=1
		exit
	end if
end do

if (noatmwfn==0) then
	write(*,"(a)") " All atom .wfn files needed have already presented in ""atomwfn"" folder, we will not calculate them"
else if (noatmwfn==1) then !Some or all atomic wfn don't exist, calc them
	!Select calculation level
	write(*,"(a)") " Note: Some or all atom .wfn files needed are not present in ""atomwfn"" folder, they must be calculated now. See Section 3.7.3 of the manual for detail."
	write(*,"(a)") " Now please input the level for calculating atom wfn files, theoretical method is optional."
	write(*,"(a)") " For example: 6-31G* or B3LYP/def2SVP    You can also add other keywords at the same time, e.g. M062X/6-311G(2d,p) scf=xqc int=ultrafine"
	read(*,"(a)") basisset !Note: 6d 10f is not required for generating wfn files, since the work has been done in L607 internally
	!Check Gaussian path
	inquire(file=gaupath,exist=alive)
	if (.not.alive) then
		write(*,*) "Could not find Gaussian path defined in ""gaupath"" variable in settings.ini"
		if (isys==1) write(*,*) "Input the path of Gaussian executable file, e.g. ""D:\study\g16w\g16.exe"""
		if (isys==2) write(*,*) "Input the path of Gaussian executable file, e.g. ""/sob/g16/g16"""
		do while(.true.)
			read(*,"(a)") gaupath
			inquire(file=gaupath,exist=alive)
			if (alive) exit
			write(*,*) "Could not find Gaussian executable file, input again"
		end do
	end if
end if

!Generate .gjf file for all elements, regardless if their wfn file have already presented, meanwhile count the total number of elements
itype=0
do i=1,nfragatm
	inquire(file=trim(tmpdir)//a(fragatm(i))%name//".gjf",exist=alive)
	if (.not.alive) then
		itype=itype+1 !How many different types
		typename(itype)=a(fragatm(i))%name
				
		if (a_org(fragatm_org(i))%index>36) then
			inquire(file=trim(tmpdir)//a(fragatm(i))%name//".wfn",exist=alive)
			if (.not.alive) then !The wfn file of the heavy element hasn't been provided in "atomwfn" and hence cannot be found in "wfntmp" here
				write(*,"(a,a,a)") " Error: Multiwfn cannot invoke Gaussian to generate wavefunction file and sphericalize density for ",a(fragatm(i))%name,", since its &
				index is larger than 36! You should provide corresponding atom .wfn files in ""atomwfn"" folder manually"
				write(*,*) "Press ENTER button to continue"
				read(*,*)
				return
			end if
		end if
		
		open(14,file=trim(tmpdir)//a(fragatm(i))%name//".gjf",status="replace")
		!If user inputted including "/" e.g. B3LYP/6-31g*, will replace default theoretical method
		if (index(basisset,'/')==0) then
			if (a(fragatm(i))%index<=20.or.a(fragatm(i))%index>=31) then
				write(14,"(a,/)") "#T out=wfn ROHF/"//trim(basisset) !Main group elements. If not use scf=sp, in g09, RO calculations for IIIA elements are to converge
				write(14,"(a,/)") "Temporary file for promolecule, ROHF"//trim(basisset)
			else
				write(14,"(a,/)") "#T out=wfn UB3LYP/"//trim(basisset) !Transition metals
				write(14,"(a,/)") "Temporary file for promolecule, UB3LYP"//trim(basisset)
			end if
		else
			if (a(fragatm(i))%index<=20.or.a(fragatm(i))%index>=31) then
				write(14,"(a,/)") "#T out=wfn RO"//trim(basisset) !Main group elements
				write(14,"(a,/)") "Temporary file for promolecule, RO"//trim(basisset)
			else
				write(14,"(a,/)") "#T out=wfn U"//trim(basisset) !Transition metals (RO may leads to convergence problem)
				write(14,"(a,/)") "Temporary file for promolecule, U"//trim(basisset)
			end if
		end if

		!Currently support up to the fourth row
		if (a(fragatm(i))%name=="H ".or.a(fragatm(i))%name=="Li".or.a(fragatm(i))%name=="Na".or.a(fragatm(i))%name=="K") write(14,*) "0 2"
		if (a(fragatm(i))%name=="Be".or.a(fragatm(i))%name=="Mg".or.a(fragatm(i))%name=="Ca") write(14,*) "0 1"
		if (a(fragatm(i))%name=="B ".or.a(fragatm(i))%name=="Al".or.a(fragatm(i))%name=="Ga") write(14,*) "0 2"
		if (a(fragatm(i))%name=="C ".or.a(fragatm(i))%name=="Si".or.a(fragatm(i))%name=="Ge") then
			if (SpherIVgroup==0) write(14,*) "0 5"
			if (SpherIVgroup==1) write(14,*) "0 3"
		end if
		if (a(fragatm(i))%name=="N ".or.a(fragatm(i))%name=="P ".or.a(fragatm(i))%name=="As") write(14,*) "0 4"
		if (a(fragatm(i))%name=="O ".or.a(fragatm(i))%name=="S ".or.a(fragatm(i))%name=="Se") write(14,*) "0 3"
		if (a(fragatm(i))%name=="F ".or.a(fragatm(i))%name=="Cl".or.a(fragatm(i))%name=="Br") write(14,*) "0 2"
		if (a(fragatm(i))%name=="He".or.a(fragatm(i))%name=="Ne".or.a(fragatm(i))%name=="Ar".or.a(fragatm(i))%name=="Kr") write(14,*) "0 1"
		if (a(fragatm(i))%name=="Sc") write(14,*) "0 2" !3d1 4s2
		if (a(fragatm(i))%name=="Ti") write(14,*) "0 3" !3d2 4s2
		if (a(fragatm(i))%name=="V ") write(14,*) "0 4" !3d3 4s2
		if (a(fragatm(i))%name=="Cr") write(14,*) "0 7" !3d5 4s1, needn't sphericalization
		if (a(fragatm(i))%name=="Mn") write(14,*) "0 6" !3d5 4s2, needn't sphericalization
		if (a(fragatm(i))%name=="Fe") write(14,*) "0 5" !3d6 4s2
		if (a(fragatm(i))%name=="Co") write(14,*) "0 4" !3d7 4s2
		if (a(fragatm(i))%name=="Ni") write(14,*) "0 3" !3d8 4s2
		if (a(fragatm(i))%name=="Cu") write(14,*) "0 2" !3d10 4s1, needn't sphericalization
		if (a(fragatm(i))%name=="Zn") write(14,*) "0 1" !3d10 4s2, needn't sphericalization
		write(14,*) a(fragatm(i))%name,0.0,0.0,0.0
		write(14,*)
		write(14,*) trim(tmpdir)//a(fragatm(i))%name//".wfn" !The output path of wfn file
		write(14,*)
		write(14,*)
		close(14)
	end if
end do

if (noatmwfn==0) then
	call delfile(trim(tmpdir)//"*.gjf") !The .gjf generated have valueless now, delete them for avoiding user's misunderstanding
else if (noatmwfn==1) then !Some wfn needs to be genereated by Gaussian and sphericalized here
	do i=1,nfragatm
		nametmp=a_org(fragatm_org(i))%name
		inquire(file=trim(tmpdir)//nametmp//".wfn",exist=alive)
		if (alive) cycle !If the .wfn file had copied from atomwfn folder, needn't recalculate

		call runcommand('"'//trim(gaupath)//'" "'//trim(tmpdir)//nametmp//'.gjf" "'//trim(tmpdir)//nametmp//'"')
		!Check if Gaussian task was successfully finished
		if (isys==1) inquire(file=trim(tmpdir)//trim(nametmp)//".out",exist=alivegauout)
		if (isys==2) inquire(file=trim(tmpdir)//trim(nametmp)//".log",exist=alivegauout)
		if (alivegauout) then
			if (isys==1) open(10,file=trim(tmpdir)//trim(nametmp)//".out",status="old")
			if (isys==2) open(10,file=trim(tmpdir)//trim(nametmp)//".log",status="old")
			call loclabel(10,"Normal termination",igaunormal)
			close(10)
			if (igaunormal==0) then
				write(*,"(a)") " Gaussian running may be failed! Please manually check Gaussian input and output files in wfntmp folder. Press ENTER button to continue"
				read(*,*)
			else if (igaunormal==1) then
				write(*,*) "Finished successfully!"
			end if
		else
			write(*,"(a)") " Gaussian running may be failed! Please manually check Gaussian input and output files in wfntmp folder"
			read(*,*)
		end if
	
		!Load and sphericalize electron density for the just generated wfn, and then save
		if (ispheratm==1.and.igaunormal==1) then
			call dealloall
			call readwfn(trim(tmpdir)//nametmp//".wfn",1)
			!Main group, restrict open-shell
			if (nametmp=="H ".or.nametmp=="Li".or.nametmp=="Na".or.nametmp=="K") MOocc(nmo)=1.0D0
			if (nametmp=="B ".or.nametmp=="Al".or.nametmp=="Ga") then
				nmo=nmo+2
				call resizebynmo(nmo,nprims) !Enlarge nmo by 2, but don't interfere nprims
				MOene(nmo-1:nmo)=MOene(nmo-2)
				MOtype(nmo-2:nmo)=1 !actually no use, because we only use atomic wfn. files to get total density
				MOocc(nmo-2:nmo)=1D0/3D0
				call orbcoeffrotate(nmo-2) !XYZ->ZXY, note: nmo-2 is original single occupied orbital
				CO(nmo-1,:)=CO(nmo-2,:)
				call orbcoeffrotate(nmo-2) !ZXY->YZX
				CO(nmo,:)=CO(nmo-2,:)
				call orbcoeffrotate(nmo-2) !YZX->XYZ, namely recovered
				!Now nmo-2,nmo-1,nmo correspond XYZ,ZXY,YZX
			end if
			if (nametmp=="C ".or.nametmp=="Si".or.nametmp=="Ge") then
				if (SpherIVgroup==0) then
					MOocc(nmo-3:nmo)=1.0D0
				else if (SpherIVgroup==1) then
					nmo=nmo+1
					call resizebynmo(nmo,nprims)
					MOene(nmo)=MOene(nmo-2) !MOene(nmo-1) is degenerate to MOene(nmo-2)
					MOtype(nmo-2:nmo)=1
					MOocc(nmo-2:nmo)=2D0/3D0
					!Rotate and copy the first occupied p orbital (nmo-5)
					call orbcoeffrotate(nmo-2) !XYZ->ZXY
					CO(nmo-1,:)=CO(nmo-2,:) !Overlap the already occupied orbital
					call orbcoeffrotate(nmo-2) !ZXY->YZX
					CO(nmo,:)=CO(nmo-2,:)
					call orbcoeffrotate(nmo-2) !YZX->XYZ, namely recovered
				end if
			end if
			if (nametmp=="N ".or.nametmp=="P ".or.nametmp=="As") MOocc(nmo-2:nmo)=1.0D0
			if (nametmp=="O ".or.nametmp=="S ".or.nametmp=="Se") MOocc(nmo-2:nmo)=4D0/3D0
			if (nametmp=="F ".or.nametmp=="Cl".or.nametmp=="Br") MOocc(nmo-2:nmo)=5D0/3D0
			!Transition metals, unrestrict open-shell, find boundary of alpha and beta first
			do ibound=2,nmo
				if (MOene(ibound)<MOene(ibound-1)) exit !from ii is beta orbitals
			end do
			!For Sc, Ti and V, rotate and duplicate d orbitals in each diection to get *near* spherical density, as for III main group
			!Note: Don't use Hartree-Fock, because correct energy sequence couldn't be reproduced, so can't be sphericalized correctly!
			if (nametmp=="Sc".or.nametmp=="Ti".or.nametmp=="V ") then !3d1 4s2, 3d2 4s2, 3d3 4s2
				if (nametmp=="Sc") then
					ibeg=1 !alpha 4s orbital, because this s orbital shows very strong unequlitity
					iend=2
				else if (nametmp=="Ti") then
					ibeg=2
					iend=3
				else if (nametmp=="V") then
					ibeg=2
					iend=4
				end if
				ienlarge=(iend-ibeg+1)*2
				call resizebynmo(nmo+ienlarge,nprims)
				ipass=0
				do iavgorb=ibeg,iend
					call orbcoeffrotate(ibound-iavgorb) !rotate this orbital
					CO(nmo+1+ipass,:)=CO(ibound-iavgorb,:) !Duplicate this orbital
					call orbcoeffrotate(ibound-iavgorb)
					CO(nmo+2+ipass,:)=CO(ibound-iavgorb,:)
					call orbcoeffrotate(ibound-iavgorb) !recover
					MOocc(ibound-iavgorb)=1D0/3D0
					MOene(nmo+1+ipass:nmo+2+ipass)=MOene(ibound-iavgorb)
					ipass=ipass+2 !next time skip nmo+1 and nmo+2
				end do
				MOocc(nmo+1:nmo+ienlarge)=1D0/3D0
				MOtype(nmo+1:nmo+ienlarge)=1
				nmo=nmo+ienlarge
			else if (nametmp=="Fe") then !3d6 4s2
				MOocc(nmo-1)=0D0 !delete the only d-beta orbital, the "nmo"th orbital is 4s-beta
				MOocc(ibound-6:ibound-2)=1.2D0 !Scatter one electron in beta-d orbital to alpha orbitals evenly. MOocc(ibound) is 4s orbital
			else if (nametmp=="Co") then !3d7 4s2
				MOocc(nmo-2:nmo-1)=0D0
				MOocc(ibound-6:ibound-2)=1.4D0
			else if (nametmp=="Ni") then !3d8 4s2
				MOocc(nmo-3:nmo-1)=0D0
				MOocc(ibound-6:ibound-2)=1.6D0
			end if
			call outwfn(trim(tmpdir)//nametmp//".wfn",0,0,10)
		end if
	end do
end if
write(*,*)

!Setup custom operation array with current size
ncustommap=nfragatm_org
if (allocated(custommapname)) deallocate(custommapname)
if (allocated(customop)) deallocate(customop)
allocate(custommapname(ncustommap))
allocate(customop(ncustommap))

!Generate atomic wfn file from element wfn file, meanwhile take them into custom operation list
do i=1,itype !Scan each atomtype in current system
	call dealloall
	call readwfn(trim(tmpdir)//typename(i)//".wfn",1)
	do j=1,nfragatm_org
		if (a_org(fragatm_org(j))%name==typename(i)) then !Find atoms attributed to current element
			a(1)%x=a_org(fragatm_org(j))%x !Modify the atomic .wfn, then output to new .wfn
			a(1)%y=a_org(fragatm_org(j))%y
			a(1)%z=a_org(fragatm_org(j))%z
			write(outwfnname,"(a2,i4,a4)") typename(i),fragatm_org(j),".wfn"
			call outwfn(trim(tmpdir)//outwfnname,0,0,10)
			custommapname(j)=trim(tmpdir)//outwfnname !Sequence is identical to atom in fragment
		end if
	end do
end do

call dealloall
call readinfile(firstfilename,1)
end subroutine



!!------- Generate a promolecular wavefunction by calculating and then combining atomic .wfn files, store to global arrays with _pmol
!  The basis set for evaluating atoms must be identical to present system  
!  Mainly used to calculate information gain at a batch of points, which needs evaluation of promolecular density frequently
!  The density calculated by this promolecular .wfn file is exactly identical to superposition of densities corresponding to all atomic .wfn files
!Since spin density is not interest in this context, spin flip is not taken into account
subroutine generate_promolwfn
use defvar
implicit real*8 (a-h,o-z)

!Generate atomic .wfn files
call setPromol

nmo_pmol=0
if (allocated(MOocc_pmol)) deallocate(MOocc_pmol,MOtype_pmol,MOene_pmol,CO_pmol)
allocate(MOocc_pmol(2*nmo),MOtype_pmol(2*nmo),MOene_pmol(2*nmo),CO_pmol(2*nmo,nprims)) !2*nmo is large enough for storing combined MOs
MOocc_pmol=0
MOene_pmol=0
MOtype_pmol=0
CO_pmol=0

!Cycling each atomic .wfn file and merge into combined wavefunction
iGTF=1
imo_pmol=1
do iatm=1,ncustommap
    call dealloall
    write(*,"(a)") " Dealing with "//trim(custommapname(iatm))
    call readinfile(custommapname(iatm),1)
    if (iGTF+nprims-1>nprims_org) then
        write(*,"(a)") " Error: The basis set used for calculating atoms must be different to that originally used for calculating molecule!"
        write(*,*) "This situation is not supported. Press ENTER button to exit"
        read(*,*)
        stop
    end if
    MOocc_pmol(imo_pmol:imo_pmol+nmo-1)=MOocc(:)
    MOene_pmol(imo_pmol:imo_pmol+nmo-1)=MOene(:)
    MOtype_pmol(imo_pmol:imo_pmol+nmo-1)=MOtype(:)
    CO_pmol(imo_pmol:imo_pmol+nmo-1,iGTF:iGTF+nprims-1)=CO(:,:)
    iGTF=iGTF+nprims
    imo_pmol=imo_pmol+nmo
end do
nmo_pmol=imo_pmol-1

write(*,*) "Done! Promolecular wavefunction has been successfully generated!"
write(*,"(a)") " Reloading "//trim(firstfilename)
call dealloall
call readinfile(firstfilename,1)

!Test if rho can be correctly calculated
!deallocate(MOocc,MOtype,MOene,CO)
!allocate(MOocc(nmo_pmol),MOene(nmo_pmol),MOtype(nmo_pmol),CO(nmo_pmol,nprims))
!nmo=nmo_pmol
!MOocc=MOocc_pmol
!MOene=MOene_pmol
!MOtype=MOtype_pmol
!CO=CO_pmol
end subroutine



!!!------------------ Generate density matrix, can be used when basis function information is available
!PS: density matrics can also be directly loaded from .fch via subroutine readfchdensmat
subroutine genP
use defvar
use util
implicit real*8 (a-h,o-z)
if (allocated(Ptot)) deallocate(Ptot)
if (allocated(Palpha)) deallocate(Palpha)
if (allocated(Pbeta)) deallocate(Pbeta)
allocate(Ptot(nbasis,nbasis))
Ptot=0
if (wfntype==1.or.wfntype==2.or.wfntype==4) then !open-shell
	allocate(Palpha(nbasis,nbasis))
	allocate(Pbeta(nbasis,nbasis))
	Palpha=0D0
	Pbeta=0D0
end if

!For SCF wavefunction, if the wavefunction has not been modified (imodwfn==0), use fast way to construct it
!However, if the wavefunction has been modified, the case may be complicated, for example, there is a hole orbital. In these cases
!We use general way (as used for post-HF) to construct density matrix
if (wfntype==0.and.imodwfn==0) then !RHF
	!Ptot=2*matmul(CObasa(:,1:nint(naelec)),transpose(CObasa(:,1:nint(naelec))))
    Ptot=matmul_blas(CObasa(:,1:nint(naelec)),CObasa(:,1:nint(naelec)),nbasis,nbasis,0,1) !Parallel MKL, fastest
    !call matprod(2,Ptot,CObasa(:,1:nint(naelec)),CObasa(:,1:nint(naelec)))
    Ptot=Ptot*2
else if (wfntype==1.and.imodwfn==0) then !UHF
	!Palpha=matmul(CObasa(:,1:nint(naelec)),transpose(CObasa(:,1:nint(naelec))))
	!Pbeta=matmul(CObasb(:,1:nint(nbelec)),transpose(CObasb(:,1:nint(nbelec))))
    Palpha=matmul_blas(CObasa(:,1:nint(naelec)),CObasa(:,1:nint(naelec)),nbasis,nbasis,0,1)
    Pbeta=matmul_blas(CObasb(:,1:nint(nbelec)),CObasb(:,1:nint(nbelec)),nbasis,nbasis,0,1)
    !call matprod(2,Palpha,CObasa(:,1:nint(naelec)),CObasa(:,1:nint(naelec)))
    !call matprod(2,Pbeta,CObasa(:,1:nint(nbelec)),CObasa(:,1:nint(nbelec)))
	Ptot=Palpha+Pbeta
else if (wfntype==2.and.imodwfn==0) then !ROHF
	Palpha=matmul(CObasa(:,1:nint(naelec)),transpose(CObasa(:,1:nint(naelec))))
	Pbeta=matmul(CObasa(:,1:nint(nbelec)),transpose(CObasa(:,1:nint(nbelec))))
	Ptot=Palpha+Pbeta
else if (wfntype==3.or.((wfntype==0.or.wfntype==2).and.imodwfn==1)) then !Restricted post-HF
	do imo=1,nmo
		if (MOocc(imo)==0D0) cycle
		Ptot=Ptot+MOocc(imo)*matmul(CObasa(:,imo:imo),transpose(CObasa(:,imo:imo)))
	end do
else if (wfntype==4.or.(wfntype==1.and.imodwfn==1)) then !Unrestricted post-HF
	do imo=1,nbasis
		if (MOocc(imo)==0D0) cycle
		Palpha=Palpha+MOocc(imo)*matmul(CObasa(:,imo:imo),transpose(CObasa(:,imo:imo)))
	end do
	do imo=1,nbasis
		if (MOocc(imo+nbasis)==0D0) cycle
		Pbeta=Pbeta+MOocc(imo+nbasis)*matmul(CObasb(:,imo:imo),transpose(CObasb(:,imo:imo)))
	end do
	Ptot=Palpha+Pbeta
end if

end subroutine



!!!------ Generate density matrix based on GTF
subroutine genPprim
use defvar
use util
implicit real*8 (a-h,o-z)

if (allocated(Ptot_prim)) deallocate(Ptot_prim)
allocate(Ptot_prim(nprims,nprims))
Ptot_prim=0

if (wfntype==1.or.wfntype==2.or.wfntype==4) then
    if (allocated(Palpha_prim)) deallocate(Palpha_prim,Pbeta_prim)
    allocate(Palpha_prim(nprims,nprims),Pbeta_prim(nprims,nprims))
    Palpha_prim=0
    Pbeta_prim=0
end if

naint=nint(naelec)
nbint=nint(nbelec)
if (wfntype==0.and.imodwfn==0) then !R wavefunction
    Ptot_prim=2*matmul_blas(transpose(CO(1:naint,:)),CO(1:naint,:),nprims,nprims,0,0)
    
else if (wfntype==1.and.imodwfn==0) then !U wavefunction
    do ibbeg=1,nmo
	    if (MOtype(ibbeg)==2) exit
    end do
    ibend=ibbeg-1+nbint
    Palpha_prim=matmul_blas(transpose(CO(1:naint,:)),CO(1:naint,:),nprims,nprims,0,0)
    if (nbint>0) then
        Pbeta_prim=matmul_blas(transpose(CO(ibbeg:ibend,:)),CO(ibbeg:ibend,:),nprims,nprims,0,0)
    end if
    Ptot_prim=Palpha_prim+Pbeta_prim
    
else if (wfntype==2.and.imodwfn==0) then !RO wavefunction
    Palpha_prim=matmul_blas(transpose(CO(1:naint,:)),CO(1:naint,:),nprims,nprims,0,0)
    if (nbint>0) then
        Pbeta_prim=matmul_blas(transpose(CO(1:nbint,:)),CO(1:nbint,:),nprims,nprims,0,0)
    end if
    Ptot_prim=Palpha_prim+Pbeta_prim
    
else if (wfntype==3.or.((wfntype==0.or.wfntype==2).and.imodwfn==1)) then !Restricted post-HF
    !$OMP PARALLEL DO SHARED(Ptot_prim) PRIVATE(iGTF,jGTF,imo) schedule(auto) NUM_THREADS(nthreads)
    do iGTF=1,nprims
        do jGTF=1,nprims
            do imo=1,nmo
                Ptot_prim(iGTF,jGTF) = Ptot_prim(iGTF,jGTF)+MOocc(imo)*CO(imo,iGTF)*CO(imo,jGTF)
            end do
        end do
    end do
    !$OMP END PARALLEL DO
    
else if (wfntype==4.or.(wfntype==1.and.imodwfn==1)) then !Unrestricted post-HF
    do ibbeg=1,nmo
	    if (MOtype(ibbeg)==2) exit
    end do
    !$OMP PARALLEL DO SHARED(Palpha_prim) PRIVATE(iGTF,jGTF,imo) schedule(auto) NUM_THREADS(nthreads)
    do iGTF=1,nprims
        do jGTF=1,nprims
            do imo=1,ibbeg-1
                Palpha_prim(iGTF,jGTF) = Palpha_prim(iGTF,jGTF)+MOocc(imo)*CO(imo,iGTF)*CO(imo,jGTF)
            end do
        end do
    end do
    !$OMP END PARALLEL DO
    !$OMP PARALLEL DO SHARED(Pbeta_prim) PRIVATE(iGTF,jGTF,imo) schedule(auto) NUM_THREADS(nthreads)
    do iGTF=1,nprims
        do jGTF=1,nprims
            do imo=ibbeg,nmo
                Pbeta_prim(iGTF,jGTF) = Pbeta_prim(iGTF,jGTF)+MOocc(imo)*CO(imo,iGTF)*CO(imo,jGTF)
            end do
        end do
    end do
    !$OMP END PARALLEL DO
    Ptot_prim=Palpha_prim+Pbeta_prim
end if
end subroutine



!!!------ Show system one-electron properties based on density matrix and integral matrix between basis functions
subroutine sys1eprop
use defvar
if (.not.allocated(Sbas)) allocate(Sbas(nbasis,nbasis))
call genSbas_curr
write(*,"(' Total number of electrons:',f16.8)") sum(Ptot*Sbas)
if (.not.allocated(Tbas)) allocate(Tbas(nbasis,nbasis))
call genTbas_curr
write(*,"(' Kinetic energy:',f18.9,' a.u.')") sum(Ptot*Tbas)
if (.not.allocated(Dbas)) allocate(Dbas(3,nbasis,nbasis))
call genDbas_curr
write(*,"(' Electric dipole moment in X/Y/Z:',3f13.7,' a.u.')") sum(Ptot*Dbas(1,:,:)),sum(Ptot*Dbas(2,:,:)),sum(Ptot*Dbas(3,:,:))
if (.not.allocated(Magbas)) allocate(Magbas(3,nbasis,nbasis))
call genMagbas_curr
write(*,"(' Magnetic dipole moment in X/Y/Z:',3f13.7,' a.u.')") sum(Ptot*Magbas(1,:,:)),sum(Ptot*Magbas(2,:,:)),sum(Ptot*Magbas(3,:,:))
if (.not.allocated(Velbas)) allocate(Velbas(3,nbasis,nbasis))
call genVelbas_curr
write(*,"(' Linear momentum in X/Y/Z:       ',3f13.7,' a.u.')") sum(Ptot*Velbas(1,:,:)),sum(Ptot*Velbas(2,:,:)),sum(Ptot*Velbas(3,:,:))
end subroutine



!!!------------- Show all properties at a point
!ifuncsel: Controls the gradient and Hessian for which function
!ifileid: Output to which file destination, of course 6=screen
subroutine showptprop(inx,iny,inz,ifuncsel,ifileid)
use util
use defvar
use function
implicit real*8(a-h,o-z)
real*8 inx,iny,inz
real*8 eigvecmat(3,3),eigval(3),elegrad(3),elehess(3,3),funcgrad(3),funchess(3,3),tmparr(3,1),tmpmat(3,3),tmpgrad1(3),tmpgrad2(3)
integer ifuncsel,ifileid
if (allocated(b)) then !If loaded file contains wavefuntion information
	call gencalchessmat(2,1,inx,iny,inz,elerho,elegrad,elehess) !Generate electron density, gradient and hessian
	write(ifileid,"(' Density of all electrons:',E18.10)") elerho
	if (ipolarpara==0) then
		tmpval=fspindens(inx,iny,inz,'s')
		write(ifileid,"(' Density of Alpha electrons:',E18.10)") (elerho+tmpval)/2D0
		write(ifileid,"(' Density of Beta electrons:',E18.10)") (elerho-tmpval)/2D0
		write(ifileid,"(' Spin density of electrons:',E18.10)") tmpval
	else if (ipolarpara==1) then
		write(ifileid,"(' Spin polarization parameter function:',E18.10)") fspindens(inx,iny,inz,'s')
	end if
	valG=lagkin(inx,iny,inz,0)
	valGx=lagkin(inx,iny,inz,1)
	valGy=lagkin(inx,iny,inz,2)
	valGz=lagkin(inx,iny,inz,3)
	write(ifileid,"(' Lagrangian kinetic energy G(r):',E18.10)") valG
	write(ifileid,"(' G(r) in X,Y,Z:',3E18.10)") valGx,valGy,valGz
	valK=Hamkin(inx,iny,inz,0)
	write(ifileid,"(' Hamiltonian kinetic energy K(r):',E18.10)") valK
! 	valKx=Hamkin(inx,iny,inz,1);valKy=Hamkin(inx,iny,inz,2);valKz=Hamkin(inx,iny,inz,3)
! 	write(ifileid,"(' K(r) in X,Y,Z:',3E18.10)") valKx,valKy,valKz
	write(ifileid,"(' Potential energy density V(r):',E18.10)") -valK-valG !When without EDF, also equals to flapl(inx,iny,inz,'t')/4.0D0-2*valG
	write(ifileid,"(' Energy density E(r) or H(r):',E18.10)") -valK
	write(ifileid,"(' Laplacian of electron density:',E18.10)") laplfac*(elehess(1,1)+elehess(2,2)+elehess(3,3))
	write(ifileid,"(' Electron localization function (ELF):',E18.10)") ELF_LOL(inx,iny,inz,"ELF")
	write(ifileid,"(' Localized orbital locator (LOL):',E18.10)") ELF_LOL(inx,iny,inz,"LOL")
	write(ifileid,"(' Local information entropy:',E18.10)") infoentro(1,inx,iny,inz)
	write(ifileid,"(' Reduced density gradient (RDG):',E18.10)") fgrad(inx,iny,inz,'r')
	write(ifileid,"(' Reduced density gradient with promolecular approximation:',E18.10)") RDGprodens(inx,iny,inz)
	write(ifileid,"(' Sign(lambda2)*rho:',E18.10)") signlambda2rho(inx,iny,inz)
	write(ifileid,"(' Sign(lambda2)*rho with promolecular approximation:',E18.10)") signlambda2rho_prodens(inx,iny,inz)
	if (pairfunctype==1) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. hole for alpha, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz)
	if (pairfunctype==2) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. hole for beta, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz)
	if (pairfunctype==4) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. fac. for alpha, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz)
	if (pairfunctype==5) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. fac. for beta, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz)
	if (pairfunctype==7) write(ifileid,"(a,3f10.5,' :',E18.10)") " Exc.-corr. dens. for alpha, ref:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz)
	if (pairfunctype==8) write(ifileid,"(a,3f10.5,' :',E18.10)") " Exc.-corr. dens. for beta, ref:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz)
	write(ifileid,"(' Source function, ref.:',3f10.5,' :',E18.10)") refx,refy,refz,srcfunc(inx,iny,inz,srcfuncmode)
	if (nmo/=0) write(ifileid,"(' Wavefunction value for orbital',i10,' :',E18.10)") iorbsel,fmo(inx,iny,inz,iorbsel)
	if (iALIEdecomp==0) then
		write(ifileid,"(' Average local ionization energy (ALIE):',E18.10)") avglocion(inx,iny,inz)
	else if (iALIEdecomp==1) then
		call avglociondecomp(ifileid,inx,iny,inz)
	end if
	write(ifileid,"(' Delta-g (under promolecular approximation):',E18.10)") delta_g_promol(inx,iny,inz)
	write(ifileid,"(' Delta-g (under Hirshfeld partition):',E18.10)") delta_g_Hirsh(inx,iny,inz)
	write(ifileid,"(' User-defined real space function:',E18.10)") userfunc(inx,iny,inz)
    if (ifPBC==0) then
	    fesptmp=nucesp(inx,iny,inz)
	    if (ifiletype==4) then
		    write(ifileid,"(' ESP from atomic charges:',E18.10)") fesptmp
	    else
		    write(ifileid,"(' ESP from nuclear charges:',E18.10)") fesptmp
	    end if
	    if (ishowptESP==1) then
		    fesptmpelec=eleesp(inx,iny,inz)
		    write(ifileid,"(' ESP from electrons:',E18.10)") fesptmpelec
		    write(ifileid,"(' Total ESP:',E18.10,' a.u. (',E14.7,' eV,',E14.7,' kcal/mol)')") &
		    fesptmpelec+fesptmp,(fesptmpelec+fesptmp)*au2eV,(fesptmpelec+fesptmp)*au2kcal
            if (iskipnuc/=0) then
		        tmpval=totespskip(inx,iny,inz,iskipnuc)
                write(ifileid,"(' Total ESP without contribution from nuclear charge of &
		        atom',i6,':',/,E18.10,' a.u. (',E15.7,' eV,',E15.7,' kcal/mol)')") iskipnuc,tmpval,tmpval*au2eV,tmpval*au2kcal
            end if
        end if
    else if (ifPBC>0) then
        write(ifileid,"(a)") " Electrostatic potential (ESP) information is not shown because it cannot be calculated for periodic systems currently"
	end if
	write(ifileid,*)
	if (ifuncsel==1) then
		write(ifileid,*) "Note: Below information are for electron density"
		funchess=elehess
		funcgrad=elegrad
	else
		if (ifuncsel==3) then
            write(ifileid,*) "Note: Below information are for Laplacian of electron density"
		else if (ifuncsel==4) then
            write(ifileid,*) "Note: Below information are for value of orbital wavefunction"
		else if (ifuncsel==9) then
            write(ifileid,*) "Note: Below information are for electron localization function"
		else if (ifuncsel==10) then
            write(ifileid,*) "Note: Below information are for localized orbital locator"
		else if (ifuncsel==12) then
            write(ifileid,*) "Note: Below information are for total ESP"
		else if (ifuncsel==100) then
            write(ifileid,*) "Note: Below information are for user-defined real space function"
		else
            write(ifileid,"(a,i4)") " Note: Below information are for real space function",ifuncsel
        end if
		call gencalchessmat(2,ifuncsel,inx,iny,inz,funcvalue,funcgrad,funchess)
	end if
	write(ifileid,*)
	write(ifileid,*) "Components of gradient in x/y/z are:"
	write(ifileid,"(3E18.10)") funcgrad(1),funcgrad(2),funcgrad(3)
	write(ifileid,"(' Norm of gradient is:',E18.10)") dsqrt(sum(funcgrad**2))
	write(ifileid,*)
	write(ifileid,*) "Components of Laplacian in x/y/z are:"
	write(ifileid,"(3E18.10)") funchess(1,1),funchess(2,2),funchess(3,3)
	write(ifileid,"(' Total:',E18.10)") funchess(1,1)+funchess(2,2)+funchess(3,3)
	write(ifileid,*)
	write(ifileid,*) "Hessian matrix:"
	write(ifileid,"(3E18.10)") funchess
	!call diagmat(funchess,eigvecmat,eigval,300,1D-12)
 	call diagsymat(funchess,eigvecmat,eigval,idiagok) !More robust
 	if (idiagok/=0) write(*,*) "Note: Diagonization of Hessian matrix failed!"
	write(ifileid,"(' Eigenvalues of Hessian:',3E18.10)") eigval(1:3)
	write(ifileid,*) "Eigenvectors(columns) of Hessian:"
	write(ifileid,"(3E18.10)") ((eigvecmat(i,j),j=1,3),i=1,3)
	write(ifileid,"(' Determinant of Hessian:',E18.10)") detmat(funchess)
	if (ifuncsel==1) then !Output ellipticity for rho
		call sort(eigval)
		eigmax=eigval(3)
		eigmed=eigval(2)
		eigmin=eigval(1)
		write(ifileid,"(a,f12.6)") " Ellipticity of electron density:",eigmin/eigmed-1
		write(ifileid,"(a,f12.6)") " eta index:",abs(eigmin)/eigmax
	end if

else !Only loaded structure, use YWT promolecule density
	if (ifiletype==4) then
		write(ifileid,"(' ESP from atomic charges:',E18.10)") nucesp(inx,iny,inz)
	else
		write(ifileid,"(' ESP from nuclear charges:',E18.10)") nucesp(inx,iny,inz)
	end if
	write(ifileid,*)
	call calchessmat_prodens(inx,iny,inz,elerho,elegrad,elehess)
	write(ifileid,"(a)") " Note: The loaded file does not contain wavefunction information, below results are evaluated from promolecule density"
	write(ifileid,"(' Density of electrons:',E18.10)") elerho
	write(ifileid,"(' Reduced density gradient:',E18.10)") RDGprodens(inx,iny,inz)
	write(ifileid,"(' Sign(lambda2)*rho:',E18.10)") signlambda2rho_prodens(inx,iny,inz)
	write(ifileid,"(' User-defined real space function:',E18.10)") userfunc(inx,iny,inz)
	write(ifileid,*)
	write(ifileid,*) "Components of gradient in x/y/z are:"
	write(ifileid,"(3E18.10)") elegrad(1),elegrad(2),elegrad(3)
	write(ifileid,"(' Norm of gradient is:',E18.10)") dsqrt(sum(elegrad**2))
	write(ifileid,*)
	write(ifileid,*) "Components of Laplacian in x/y/z are:"
	write(ifileid,"(3E18.10)") elehess(1,1),elehess(2,2),elehess(3,3)
	write(ifileid,"(' Total:',E18.10)") elehess(1,1)+elehess(2,2)+elehess(3,3)
	write(ifileid,*)
	write(ifileid,*) "Hessian matrix:"
	write(ifileid,"(3E18.10)") elehess
	call diagmat(elehess,eigvecmat,eigval,300,1D-12)
	write(ifileid,"(' Eigenvalues of Hessian:',3E18.10)") eigval(1:3)
	write(ifileid,*) "Eigenvectors (columns) of Hessian:"
	write(ifileid,"(3E18.10)") ((eigvecmat(i,j),j=1,3),i=1,3)
end if
end subroutine



!!!------- Decompose property at a point as contribution from various orbitals
subroutine decompptprop(x,y,z)
use defvar
use util
use function
implicit real*8(a-h,o-z)
character c2000tmp*2000
real*8 MOocctmp(nmo)
real*8,allocatable :: valarr(:)
integer,allocatable :: orbidx(:)
write(*,*) "Select the function to be studied"
call funclist
read(*,*) ifunc
write(*,*) "Input range of orbitals, e.g. 3,6-8,12-15"
write(*,"(a)") " Note: If press ENTER button directly, all occupied orbitals will be taken into account, &
and 10 orbitals having largest contributions will be shown"
read(*,"(a)") c2000tmp
if (c2000tmp==" ") then
	norb=count(MOocc(1:nmo)/=0D0)
	allocate(orbidx(norb),valarr(norb))
	idx=0
	do imo=1,nmo
		if (MOocc(imo)/=0) then
			idx=idx+1
			orbidx(idx)=imo
		end if
	end do
else
	call str2arr(c2000tmp,norb)
	allocate(orbidx(norb),valarr(norb))
	call str2arr(c2000tmp,norb,orbidx)
end if

totval=calcfuncall(ifunc,x,y,z)
MOocctmp=MOocc
sumval=0
do itmp=1,norb
	iorb=orbidx(itmp)
	MOocc=0
	MOocc(iorb)=MOocctmp(iorb)
	valarr(itmp)=calcfuncall(ifunc,x,y,z)
	sumval=sumval+valarr(itmp)
end do
MOocc=MOocctmp

call sortr8(valarr,"abs",orbidx)
call invarrr8(valarr,orbidx)

nout=norb
if (c2000tmp==" ".and.norb>10) nout=10
do itmp=1,nout
	iorb=orbidx(itmp)
	if (ifunc==1) then !For electron density, also print percentage contribution
		write(*,"(' Contribution from orbital',i6,' (occ=',f9.6,'):',f14.6,' a.u. (',f6.2,'% )')") iorb,MOocc(iorb),valarr(itmp),valarr(itmp)/totval*100
	else
		write(*,"(' Contribution from orbital',i6,' (occ=',f9.6,'):',1E16.8,' a.u.')") iorb,MOocc(iorb),valarr(itmp)
	end if
end do

if (ifunc==1) then
	write(*,"(' Sum of above values:',f16.8,' a.u. ( ',f6.2,'% )')") sumval,sumval/totval*100
	write(*,"(' Exact value:',f16.8,' a.u.')") totval
else
	write(*,"(' Sum of above values:',1E16.8,' a.u.')") sumval
	write(*,"(' Exact value:',1E16.8,' a.u.')") totval
end if
end subroutine



!!!------------------- Delete virtual orbitals higher than LUMO+10 for HF/DFT wavefunctions
!Each time "delvirob" has been called, "delvirorb_back" should be then called to recover previous wavefunction status
!infomode=1 means show prompt
subroutine delvirorb(infomode)
use defvar
implicit real*8 (a-h,o-z)
integer :: infomode,nvirsave=10 !Lowest "nvirsave" virtual orbitals will be reserved
if (.not.allocated(CObasa)) return !Only works when all orbitals are available, which implies CObasa is allocated
if (idelvirorb==0) return !Do not delete orbitals as explicitly requested in settings.ini
!Linear response kernel, local electron affinity, orbital-weighted Fukui/dual descriptor require all orbital information
if (iuserfunc==24.or.iuserfunc==27.or.iuserfunc==28.or.iuserfunc==29.or.iuserfunc==95.or.iuserfunc==96.or.iuserfunc==97.or.iuserfunc==98) return
!if (imodwfn==1) return !Wavefunction has been modified, don't make thing more complicated!
if (wfntype==3.or.wfntype==4) return !This routine doesn't work for post-HF cases
if (ifdelvirorb==1) return !This routine has already been called while delvirorb_back was not used after that

if (allocated(CO_back)) deallocate(CO_back,MOene_back,MOocc_back,MOtype_back)
allocate(CO_back(nmo,nprims),MOene_back(nmo),MOocc_back(nmo),MOtype_back(nmo))
nmo_back=nmo
CO_back=CO
MOene_back=MOene
MOocc_back=MOocc
MOtype_back=MOtype

call getHOMOidx !Return idxHOMO (R, RO), or idxHOMO and idxHOMOb (U)
if (wfntype==0.or.wfntype==2) then !RHF, ROHF
	if (nmo<=idxHOMO+nvirsave) return
	nmo=idxHOMO+nvirsave !Simply shield those virtual orbitals
else if (wfntype==1) then !Reserve up to LUMO+10 for alpha, and same number of orbitals for beta
	if (nmo/2<=idxHOMO+nvirsave) return !naelec is always >= nbelec
	nperserve=idxHOMO+nvirsave !Number of lowest preserved alpha orbitals
	!CObasa and CObasb needn't to be modified, because they are not directly involved in real space function calculation
	CO(nperserve+1:2*nperserve,:)=CO(nmo/2+1:nmo/2+nperserve,:)
	MOene(nperserve+1:2*nperserve)=MOene(nmo/2+1:nmo/2+nperserve)
	MOocc(nperserve+1:2*nperserve)=MOocc(nmo/2+1:nmo/2+nperserve)
	MOtype(nperserve+1:2*nperserve)=MOtype(nmo/2+1:nmo/2+nperserve)
	nmo=2*nperserve
end if

ifdelvirorb=1 !delvirorb has taken effect
if (infomode==1) write(*,"(a,/)") " Note: Virtual orbitals higher than LUMO+9 have been temporarily discarded for saving computational time"
end subroutine



!!!------------------- Recover status prior to calling delvirorb
!infomode=1 means show prompt
subroutine delvirorb_back(infomode)
use defvar
implicit real*8 (a-h,o-z)
integer infomode
if (ifdelvirorb==1) then
    nmo=nmo_back
    CO=CO_back
    MOene=MOene_back
    MOocc=MOocc_back
    MOtype=MOtype_back
    deallocate(CO_back,MOene_back,MOocc_back,MOtype_back)
    ifdelvirorb=0 !Has been restored
    if (infomode==1) write(*,"(a)") " Note: Previous orbital information has been restored"
end if
end subroutine



!!----------- Generate unique GTF information (b_uniq, CO_uniq)
!This reduce number of GTFs if generally contracted basis set is used, then cost in subroutine "orbserv" can be lowered
!infomode=0: Print number of unique GTFs. =1: Do not print
!Note that after using it, del_GTFuniq should be invoked as early as possible, otherwise when CO is changed due to some reasons, CO_uniq will be out-of-dated
subroutine gen_GTFuniq(infomode)
use defvar
integer ifcombined(nprims)

if (.not.allocated(b)) return
if (allocated(b_uniq)) deallocate(b_uniq)
if (allocated(CO_uniq)) deallocate(CO_uniq)

do itime=1,2
	nprims_uniq=0
	ifcombined=0
	do iGTF=1,nprims
		if (ifcombined(iGTF)==0) then
			nprims_uniq=nprims_uniq+1
            if (itime==2) then
				b_uniq(nprims_uniq)=b(iGTF)
				CO_uniq(:,nprims_uniq)=CO(:,iGTF)
            end if
			do jGTF=iGTF+1,nprims
				if (b(jGTF)%center==b(iGTF)%center.and.b(jGTF)%type==b(iGTF)%type.and.b(jGTF)%exp==b(iGTF)%exp) then
					if (itime==2) CO_uniq(:,nprims_uniq)=CO_uniq(:,nprims_uniq)+CO(:,jGTF)
					ifcombined(jGTF)=1
				end if
			end do
        end if
	end do
    if (itime==1) allocate(b_uniq(nprims_uniq),CO_uniq(nmo,nprims_uniq))
end do

if (infomode==0) write(*,"(' Unique GTFs have been constructed. Number of unique GTFs:',i6)") nprims_uniq
end subroutine



!!----------- Destory unique GTF information (b_uniq, CO_uniq)
subroutine del_GTFuniq
use defvar
nprims_uniq=0
if (allocated(b_uniq)) deallocate(b_uniq)
if (allocated(CO_uniq)) deallocate(CO_uniq)
end subroutine



!!!-------- imode=1: Convert unit of grid/plane parameters from Bohr to Angstrom. =2: Convert them back
subroutine convgridlenunit(imode)
use defvar
implicit none
integer imode
real*8 scll
if (imode==1) scll=b2a
if (imode==2) scll=1/b2a
orgx=orgx*scll
orgy=orgy*scll
orgz=orgz*scll
orgx2D=orgx2D*scll
orgy2D=orgy2D*scll
orgz2D=orgz2D*scll
endx=endx*scll
endy=endy*scll
endz=endz*scll
dx=dx*scll;gridv1=gridv1*scll
dy=dy*scll;gridv2=gridv2*scll
dz=dz*scll;gridv3=gridv3*scll
v1x=v1x*scll
v1y=v1y*scll
v1z=v1z*scll
v2x=v2x*scll
v2y=v2y*scll
v2z=v2z*scll
a1x=a1x*scll
a1y=a1y*scll
a1z=a1z*scll
a2x=a2x*scll
a2y=a2y*scll
a2z=a2z*scll
a3x=a3x*scll
a3y=a3y*scll
a3z=a3z*scll
d1=d1*scll
d2=d2*scll
end subroutine



!!-------- Deallocate all arrays about wavefunction except for the _org ones, and re-initialize some variables
subroutine dealloall
use defvar

call delvirorb_back(0) !If delvirorb has taken effect, use this routine to deallocate relevant arrays
call del_GTFuniq
if (allocated(a)) deallocate(a)
if (allocated(b)) deallocate(b)
if (allocated(CO)) deallocate(CO)
if (allocated(MOocc)) deallocate(MOocc)
if (allocated(MOsym)) deallocate(MOsym)
if (allocated(MOene)) deallocate(MOene)
if (allocated(MOtype)) deallocate(MOtype)
if (allocated(b_EDF)) then
	deallocate(CO_EDF,b_EDF)
	nEDFprims=0
	nEDFelec=0
end if
if (allocated(connmat)) deallocate(connmat)
!Related to basis functions
if (allocated(shtype)) deallocate(shtype,shcen,shcon,primshexp,primshcoeff)
if (allocated(basshell)) deallocate(basshell,bascen,bastype,basstart,basend,primstart,primend,primconnorm)
if (allocated(shtypeCar)) deallocate(shtypeCar)
if (allocated(CObasa)) deallocate(CObasa)
if (allocated(CObasb)) deallocate(CObasb)
if (allocated(Ptot)) deallocate(Ptot)
if (allocated(Palpha)) deallocate(Palpha)
if (allocated(Pbeta)) deallocate(Pbeta)
if (allocated(Sbas)) deallocate(Sbas)
if (allocated(Dbas)) deallocate(Dbas)
if (allocated(DorbA)) deallocate(DorbA)
if (allocated(DorbB)) deallocate(DorbB)
if (allocated(Magbas)) deallocate(Magbas)
if (allocated(Quadbas)) deallocate(Quadbas)
if (allocated(Octobas)) deallocate(Octobas)

if (allocated(Ptot_prim)) deallocate(Ptot_prim)
if (allocated(Palpha_prim)) deallocate(Palpha_prim)
if (allocated(Pbeta_prim)) deallocate(Pbeta_prim)
if (allocated(Dprim)) deallocate(Dprim)
if (allocated(Quadprim)) deallocate(Quadprim)
if (allocated(Octoprim)) deallocate(Octoprim)

if (allocated(frag1)) deallocate(frag1)
if (allocated(frag2)) deallocate(frag2)
if (allocated(fragatm)) deallocate(fragatm)

!Arrays used by delvirorb for backing up original MO information
if (allocated(CO_back)) deallocate(CO_back,MOene_back,MOocc_back,MOtype_back)

loadmulti=-99
loadcharge=-99
totenergy=0
virialratio=2
nelec=0
naelec=0
nbelec=0
nmo=0
nprims=0
ncenter=0
iresinfo=0
nbasis=0
nindbasis=0

ifPBC=0
cellv1=0;cellv2=0;cellv3=0
if (allocated(Sbas_PBC)) deallocate(Sbas_PBC)

call deallo_excitinfo !Deallocate excited state information
call deallo_basinana(0) !Deallocate basin analysis information

!Make LIBRETA in uninitialized status
if_initlibreta=0
nstates=0
end subroutine


!!-------- Deallocate all arrays about the system for the _org ones
subroutine dealloall_org
use defvar
firstfilename=" "
ncenter_org=0
nmo_org=0
nprims_org=0
if (allocated(a_org)) deallocate(a_org)
if (allocated(b_org)) deallocate(b_org,CO_org,MOocc_org,MOene_org)
if (allocated(Sbas_org)) deallocate(Sbas_org)
if (allocated(CObasa_org)) deallocate(CObasa_org)
if (allocated(CObasb_org)) deallocate(CObasb_org)
if (allocated(fragatm_org)) deallocate(fragatm_org)
ifPBC_org=0
cellv1_org=0
cellv2_org=0
cellv3_org=0
end subroutine




!!------- Generate atomic/fragmental Hirshfeld weight and store it to planemat, calculate free-atom/fragmental density and store it to planemattmp
!The atoms in the fragment is inputted as "selatm" array, nselatm is the number of its elements
!if itype=1, use atomic wavefunction to calculate Hirshfeld weight, and setpromol must have been invoked; if =2, use built-in atomic density to generate it
subroutine genHirshplanewei(selatm,nselatm,itype)
use defvar
use function
implicit real*8 (a-h,o-z)
integer selatm(nselatm),nselatm,itype
if (allocated(planemat)) deallocate(planemat)
if (allocated(planemattmp)) deallocate(planemattmp)
allocate(planemat(ngridnum1,ngridnum2),planemattmp(ngridnum1,ngridnum2))
planemat=0D0
planemattmp=0D0
do iatm=1,ncenter_org !Calc free atomic density of each atom, get promolecular density and Hirshfeld weight of present atom
	iyes=0
	if (any(selatm==iatm)) iyes=1
	if (itype==1) then
		call dealloall
		call readwfn(custommapname(iatm),1)
	end if
	!$OMP PARALLEL DO private(i,j,rnowx,rnowy,rnowz,tmpval) shared(planemat) schedule(dynamic) NUM_THREADS(nthreads)
	do i=1,ngridnum1 !First calculate promolecular density and store it to planemat
		do j=1,ngridnum2
			rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x
			rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y
			rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z
			if (itype==1) then
				tmpval=fdens(rnowx,rnowy,rnowz)
			else
				tmpval=calcatmdens(iatm,rnowx,rnowy,rnowz,0)
			end if
			planemat(i,j)=planemat(i,j)+tmpval
			if (iyes==1) planemattmp(i,j)=planemattmp(i,j)+tmpval
		end do
	end do
	!$OMP END PARALLEL DO
end do
if (itype==1) then
	call dealloall
	call readinfile(firstfilename,1) !Retrieve the first loaded file(whole molecule)
end if

do i=1,ngridnum1 !Calculate Hirshfeld weighting function
	do j=1,ngridnum2
		if (planemat(i,j)/=0D0) then
			planemat(i,j)=planemattmp(i,j)/planemat(i,j)
		else
			planemat(i,j)=0D0
		end if
	end do
end do
end subroutine



!!----- Generate atomic Hirshfeld weight and store it to cubmat
!The atoms in the fragment is inputted as "selatm" array, nselatm is the number of its elements
!if itype=1, use atomic wavefunction to calculate Hirshfeld weight, and setpromol must have been invoked; if =2, use built-in atomic density to generate it
subroutine genHirshcubewei(selatm,nselatm,itype)
use defvar
use function
implicit real*8 (a-h,o-z)
integer selatm(nselatm),nselatm,itype
if (allocated(cubmat)) deallocate(cubmat)
if (allocated(cubmattmp)) deallocate(cubmattmp)
allocate(cubmat(nx,ny,nz),cubmattmp(nx,ny,nz))
cubmat=0D0
cubmattmp=0D0
do iatm=1,ncenter_org
	write(*,"(' Finished',i6,'  /',i6)") iatm,ncenter_org
	if (itype==1) then
		call dealloall
		call readwfn(custommapname(iatm),1)
	end if
	!$OMP PARALLEL DO SHARED(cubmat,cubmattmp,ifinish) PRIVATE(i,j,k,tmpx,tmpy,tmpz,tmpval) schedule(dynamic) NUM_THREADS(nthreads)
	do k=1,nz !First calculate promolecular density and store it to cubmat
		do j=1,ny
			do i=1,nx
                call getgridxyz(i,j,k,tmpx,tmpy,tmpz)
				if (itype==1) then
					tmpval=fdens(tmpx,tmpy,tmpz)
				else
					tmpval=calcatmdens(iatm,tmpx,tmpy,tmpz,0)
				end if
				cubmat(i,j,k)=cubmat(i,j,k)+tmpval
				if (any(selatm==iatm)) cubmattmp(i,j,k)=cubmattmp(i,j,k)+tmpval
			end do
		end do
	end do
	!$OMP END PARALLEL DO
end do
if (itype==1) then
	call dealloall
	call readinfile(firstfilename,1) !Retrieve the first loaded file(whole molecule)
end if

do k=1,nz !Calculate Hirshfeld weighting function
	do j=1,ny
		do i=1,nx
			if (cubmat(i,j,k)/=0D0) then
				cubmat(i,j,k)=cubmattmp(i,j,k)/cubmat(i,j,k)
			else
				cubmat(i,j,k)=0D0
			end if
		end do
	end do
end do
end subroutine



!!--- Generate single-center integration grid for Becke's integration. Not adapted according to element. Return iradcut and gridatm
subroutine gen1cintgrid(gridatm,iradcut)
use defvar
implicit real*8 (a-h,o-z)
integer iradcut
real*8 potx(sphpot),poty(sphpot),potz(sphpot),potw(sphpot)
type(content) gridatm(radpot*sphpot)
call Lebedevgen(sphpot,potx,poty,potz,potw)
iradcut=0 !Before where the radial points will be cut
parm=1D0
do i=1,radpot !Combine spherical point&weights with second kind Gauss-Chebyshev method for radial part
	radx=cos(i*pi/(radpot+1))
	radr=(1+radx)/(1-radx)*parm !Becke transform
	radw=2*pi/(radpot+1)*parm**3 *(1+radx)**2.5D0/(1-radx)**3.5D0 *4*pi
	gridatm( (i-1)*sphpot+1:i*sphpot )%x=radr*potx
	gridatm( (i-1)*sphpot+1:i*sphpot )%y=radr*poty
	gridatm( (i-1)*sphpot+1:i*sphpot )%z=radr*potz
	gridatm( (i-1)*sphpot+1:i*sphpot )%value=radw*potw
	if (radcut/=0D0.and.iradcut==0.and.radr<radcut) iradcut=i-1
end do
end subroutine
!!--- Generate Becke weight for a batch of points around iatm, sharpness parameter=3
!!--- Input: iatm, iradcut, gridatm   Return: beckeweigrid
subroutine gen1cbeckewei(iatm,iradcut,gridatm,beckeweigrid)
use defvar
implicit real*8 (a-h,o-z)
integer iatm,iradcut
real*8 beckeweigrid(radpot*sphpot),smat(ncenter,ncenter),Pvec(ncenter)
type(content) gridatm(radpot*sphpot)
!$OMP parallel do shared(beckeweigrid) private(i,rnowx,rnowy,rnowz,smat,ii,ri,jj,rj,rmiu,chi,uij,aij,tmps,Pvec) num_threads(nthreads) schedule(dynamic)
do i=1+iradcut*sphpot,radpot*sphpot
	smat=1D0
	rnowx=gridatm(i)%x
	rnowy=gridatm(i)%y
	rnowz=gridatm(i)%z
	do ii=1,ncenter
		ri=dsqrt( (rnowx-a(ii)%x)**2+(rnowy-a(ii)%y)**2+(rnowz-a(ii)%z)**2 )
		do jj=1,ncenter
			if (ii==jj) cycle
            !if (distmat(ii,jj)>12) cycle !This is found to be quite safe, the error is negligible, and the cost can be modestly reduced
			rj=dsqrt( (rnowx-a(jj)%x)**2+(rnowy-a(jj)%y)**2+(rnowz-a(jj)%z)**2 )
			rmiu=(ri-rj)/distmat(ii,jj)
 			!Adjust for heteronuclear
			chi=covr_tianlu(a(ii)%index)/covr_tianlu(a(jj)%index)
			uij=(chi-1)/(chi+1)
			aij=uij/(uij**2-1)
			if (aij>0.5D0) aij=0.5D0
			if (aij<-0.5D0) aij=-0.5D0
			rmiu=rmiu+aij*(1-rmiu**2)
			tmps=rmiu
			do iter=1,3
				tmps=1.5D0*(tmps)-0.5D0*(tmps)**3
			end do
			smat(ii,jj)=0.5D0*(1-tmps)
		end do
	end do
	Pvec=1D0
	do ii=1,ncenter
		Pvec=Pvec*smat(:,ii)
	end do
	beckeweigrid(i)=Pvec(iatm)/sum(Pvec)
    !write(12,"(2i6,4f12.6)") iatm,i,rnowx,rnowy,rnowz,beckeweigrid(i)
end do
!$OMP end parallel do
end subroutine



!!--------- A standalone routine to calculate atomic contribution to specific real space function. Employed by IFCT analysis
!iparttype=1: Becke partition
!atmcontri: Returned array containing atomic contribution
!ifunc: The function to be evaluated
subroutine atmcontrifunc(iparttype,atmcontri,ifunc)
use defvar
use function
implicit real*8 (a-h,o-z)
integer ifunc
type(content) gridatm(radpot*sphpot),gridatmorg(radpot*sphpot)
real*8 atmcontri(ncenter),beckeweigrid(radpot*sphpot)

atmcontri=0
! write(*,"(' Radial points:',i5,'    Angular points:',i5,'   Total:',i10,' per center')") radpot,sphpot,radpot*sphpot
if (iparttype==1) then !Becke partition
	call gen1cintgrid(gridatmorg,iradcut)
	!$OMP PARALLEL DO SHARED(atmcontri) PRIVATE(iatm,ipt,gridatm,beckeweigrid,funcval) schedule(dynamic) NUM_THREADS(nthreads)
	do iatm=1,ncenter
! 		write(*,"(' Processing center',i6,'(',a2,')   /',i6)") iatm,a(iatm)%name,ncenter
		gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule
		gridatm%y=gridatmorg%y+a(iatm)%y
		gridatm%z=gridatmorg%z+a(iatm)%z
		call gen1cbeckewei(iatm,iradcut,gridatm,beckeweigrid)
		do ipt=1+iradcut*sphpot,radpot*sphpot
			funcval=calcfuncall(ifunc,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z)
			atmcontri(iatm)=atmcontri(iatm)+funcval*gridatmorg(ipt)%value*beckeweigrid(ipt)
		end do
	end do
	!$OMP END PARALLEL DO
end if
end subroutine



!!--------- Output all Becke's integration points to intpt.txt in current folder
subroutine outBeckeintpt
use defvar
implicit real*8 (a-h,o-z)
real*8 beckeweigrid(radpot*sphpot)
type(content) gridatm(radpot*sphpot),gridatmorg(radpot*sphpot)
character outwei

write(*,*) "See Appendix 6 of Multiwfn manual for detail of this function"
write(*,*)
write(*,*) "Also export Becke's integration weights? (y/n)"
read(*,*) outwei

open(10,file="intpt.txt",status="replace")
write(*,"(' Radial points:',i5,'    Angular points:',i5,'   Total:',i10,' per center')") radpot,sphpot,radpot*sphpot
if (outwei=='y') then
    write(*,*) "Calculating and exporting, please wait..."
else
    write(*,*) "Exporting, please wait..."
end if
call gen1cintgrid(gridatmorg,iradcut)
do iatm=1,ncenter
	gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule
	gridatm%y=gridatmorg%y+a(iatm)%y
	gridatm%z=gridatmorg%z+a(iatm)%z
    if (outwei=='n') then
        do ipt=1+iradcut*sphpot,radpot*sphpot
            write(10,"(i6,4E16.8)") iatm,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z,gridatmorg(ipt)%value
        end do
    else
	    call gen1cbeckewei(iatm,iradcut,gridatm,beckeweigrid)
        do ipt=1+iradcut*sphpot,radpot*sphpot
            write(10,"(i6,5E16.8)") iatm,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z,gridatmorg(ipt)%value,beckeweigrid(ipt)
        end do
    end if
end do  
close(10)
write(*,"(a)") " Done! All Becke's integration points have been exported to intpt.txt in current folder"
write(*,*) "Column 1: Index of the atom that the point belongs to"
write(*,*) "Column 2/3/4: X/Y/Z coordinate of the point in Bohr"
write(*,"(a)") " Column 5: Single center integration weight (second kind Gauss-Chebyshev for radial part and Lebedev for angular part"
if (outwei=='y') write(*,*) "Column 6: Becke's integration weight"
end subroutine



!!--------- Convert 1RDM in MO basis outputted by MRCC program to natural orbitals
!In CCDENSITIES, the density matrix is represented in MO basis
!When frozen core is enabled, the indices are counted from the first correlated orbital
subroutine MRCC_gennatorb
use defvar
use util
character c200tmp*200,selectyn,outname*200
real*8,allocatable :: eigvecmat(:,:),eigvalarr(:),tmparr(:)
do while(.true.)
	write(*,*) "Input the path of CCDENSITIES file, e.g. C:\lovelive\CCDENSITIES"
    write(*,*) "If press ENTER button directly, CCDENSITIES in current folder will be loaded"
	read(*,"(a)") c200tmp
    if (c200tmp==" ") c200tmp="CCDENSITIES"
	inquire(file=c200tmp,exist=alive)
	if (alive) exit
	write(*,*) "Error: Cannot find the file, input again"
end do
write(*,*)
write(*,*) "Input the number of frozen orbitals, e.g. 3"
write(*,*) "If no orbitals are frozen, simply input 0"
write(*,"(a)") " PS: For unrestricted reference, if you input n, the n lowest alpha and n lowest beta MOs will be regarded as frozen"
read(*,*) nfrz
write(*,*) "Please wait..."
if (wfntype==0) then !RHF reference
	open(10,file=c200tmp,status="old")
	Ptot=0
	do while(.true.)
		read(10,*,iostat=ierror) tmp,i,j,k,l
		if (ierror/=0) exit
		if (k==0.and.l==0) then !Only load 1RDM
			Ptot(i+nfrz,j+nfrz)=tmp
			Ptot(j+nfrz,i+nfrz)=tmp
		end if
	end do
	close(10)
	do ifrz=1,nfrz
		Ptot(ifrz,ifrz)=2D0
	end do
	allocate(eigvecmat(nbasis,nbasis),eigvalarr(nbasis),tmparr(nbasis))
	call diagsymat(Ptot,eigvecmat,eigvalarr,istat)
	MOocc=eigvalarr
	!Currently the occupation is from low to high, now invert the sequence
	do i=1,int(nmo/2D0)
		idx=i
		jdx=nmo+1-i
		tmp=MOocc(idx)
		MOocc(idx)=MOocc(jdx)
		MOocc(jdx)=tmp
		tmparr=eigvecmat(:,idx)
		eigvecmat(:,idx)=eigvecmat(:,jdx)
		eigvecmat(:,jdx)=tmparr
	end do
	CObasa=matmul(CObasa,eigvecmat)
	wfntype=3
	write(*,*) "Occupation numbers:"
	write(*,"(6f12.8)") MOocc
else if (wfntype==1) then !UHF reference
	!In CCDENSITIES, the sequence is:
	!2RDM-alpha
	!  0.00000000000000000000E+00   0   0   0   0
	!2RDM-beta
	!  0.00000000000000000000E+00   0   0   0   0
	!Unknown
	!  0.00000000000000000000E+00   0   0   0   0
	!1RDM-alpha
	!  0.00000000000000000000E+00   0   0   0   0
	!1RDM-beta
	!  0.00000000000000000000E+00   0   0   0   0
	open(10,file=c200tmp,status="old")
	Palpha=0
	Pbeta=0
	itime=0
	do while(.true.)
		read(10,*) tmp,i,j,k,l
		if (i==0.and.j==0.and.k==0.and.l==0) then
			itime=itime+1
			if (itime==5) exit
			cycle
		end if
		if (itime==3) then
			Palpha(i+nfrz,j+nfrz)=tmp
			Palpha(j+nfrz,i+nfrz)=tmp
		else if (itime==4) then
			Pbeta(i+nfrz,j+nfrz)=tmp
			Pbeta(j+nfrz,i+nfrz)=tmp
		end if
	end do
	close(10)
	do ifrz=1,nfrz
		Palpha(ifrz,ifrz)=1D0
		Pbeta(ifrz,ifrz)=1D0
	end do
	allocate(eigvecmat(nbasis,nbasis),eigvalarr(nbasis),tmparr(nbasis))
	!Alpha part
	call diagsymat(Palpha,eigvecmat,eigvalarr,istat)
	MOocc(1:nbasis)=eigvalarr
	do i=1,int(nbasis/2D0)
		idx=i
		jdx=nbasis+1-i
		tmp=MOocc(idx)
		MOocc(idx)=MOocc(jdx)
		MOocc(jdx)=tmp
		tmparr=eigvecmat(:,idx)
		eigvecmat(:,idx)=eigvecmat(:,jdx)
		eigvecmat(:,jdx)=tmparr
	end do
	CObasa=matmul(CObasa,eigvecmat)
	write(*,*) "Occupation number of Alpha part:"
	write(*,"(6f12.6)") MOocc(1:nbasis)
	!Beta part
	call diagsymat(Pbeta,eigvecmat,eigvalarr,istat)
	MOocc(nbasis+1:nmo)=eigvalarr
	do i=1,int(nbasis/2D0)
		idx=nbasis+i
		jdx=nmo+1-i
		tmp=MOocc(idx)
		MOocc(idx)=MOocc(jdx)
		MOocc(jdx)=tmp
		tmparr=eigvecmat(:,i)
		eigvecmat(:,i)=eigvecmat(:,nbasis+1-i)
		eigvecmat(:,nbasis+1-i)=tmparr
	end do
	CObasb=matmul(CObasb,eigvecmat)
	write(*,*) "Occupation number of Beta part:"
	write(*,"(6f12.6)") MOocc(nbasis+1:nmo)
	wfntype=4
end if

call genP
MOene=0

call path2filename(filename,c200tmp)
outname=trim(c200tmp)//".mwfn"
call outmwfn(outname,10,0)
write(*,"(/,a)") " All natural orbitals has been exported to "//trim(outname)//" in current folder"
write(*,"(a)") " Do you want to load it now so that you can perform wavefunction analysis based on the natural orbitals at the corresponding level? (y/n)"
read(*,*) selectyn
call dealloall
if (selectyn=='y'.or.selectyn=='Y') then
    write(*,*) "Loading "//trim(outname)
    call readinfile(outname,1)
    write(*,"(a)") " Loading finished!"
else
    write(*,"(a)") " Reloading "//trim(filename)//" to recover initial status..."
    call readinfile(filename,1)
    write(*,*) "Loading finished!"
end if

end subroutine




!!----------- Generate spherical harmonic -> Cartesian basis function conversion table for d,f,g,h.
!iprog=1: for readfch;  iprog=2: for readmolden
!The table comes from IJQC,54,83, which is used by Gaussian
!The sequence of d and f shell is also identical to .molden convention, but for g, another conversion table is used, &
!since in Multiwfn g Cartesian shell starts from ZZZZ, but that of .molden starts from xxxx
subroutine gensphcartab(iprog,matd,matf,matg,math)
real*8 matd(6,5),matf(10,7),matg(15,9),math(21,11)
integer iprog
matd=0D0
matf=0D0
matg=0D0
math=0D0
! From 5D: D 0,D+1,D-1,D+2,D-2
! To 6D:  1  2  3  4  5  6
!        XX,YY,ZZ,XY,XZ,YZ
!
! D0=-0.5*XX-0.5*YY+ZZ
matd(1:3,1)=(/ -0.5D0,-0.5D0,1D0 /)
! D+1=XZ
matd(5,2)=1D0
! D-1=YZ
matd(6,3)=1D0
! D+2=SQRT(3)/2*(XX-YY)
matd(1:2,4)=(/ sqrt(3D0)/2D0,-sqrt(3D0)/2D0 /)
! D-2=XY
matd(4,5)=1D0

! From 7F: F 0,F+1,F-1,F+2,F-2,F+3,F-3
! To 10F:  1   2   3   4   5   6   7   8   9  10      
!         XXX,YYY,ZZZ,XYY,XXY,XXZ,XZZ,YZZ,YYZ,XYZ (Gaussian sequence, not identical to Multiwfn)
!
! F 0=-3/(2*√5)*(XXZ+YYZ)+ZZZ
matf(3,1)=1D0
matf(6,1)=-1.5D0/sqrt(5D0)
matf(9,1)=-1.5D0/sqrt(5D0)
! F+1=-√(3/8)*XXX-√(3/40)*XYY+√(6/5)*XZZ
matf(1,2)=-sqrt(3D0/8D0)
matf(4,2)=-sqrt(3D0/40D0)
matf(7,2)=sqrt(6D0/5D0)
! F-1=-√(3/40)*XXY-√(3/8)*YYY+√(6/5)*YZZ
matf(2,3)=-sqrt(3D0/8D0)
matf(5,3)=-sqrt(3D0/40D0)
matf(8,3)=sqrt(6D0/5D0)
! F+2=√3/2*(XXZ-YYZ)
matf(6,4)=sqrt(3D0)/2D0
matf(9,4)=-sqrt(3D0)/2D0
! F-2=XYZ
matf(10,5)=1D0
! F+3=√(5/8)*XXX-3/√8*XYY
matf(1,6)=sqrt(5D0/8D0)
matf(4,6)=-3D0/sqrt(8D0)
! F-3=3/√8*XXY-√(5/8)*YYY
matf(2,7)=-sqrt(5D0/8D0)
matf(5,7)=3D0/sqrt(8D0)

if (iprog==1) then !for .fch
	! From 9G: G 0,G+1,G-1,G+2,G-2,G+3,G-3,G+4,G-4
	! To 15G:   1    2    3    4    5    6    7    8
	!         ZZZZ,YZZZ,YYZZ,YYYZ,YYYY,XZZZ,XYZZ,XYYZ
	!           9   10   11   12   13   14   15
	!         XYYY,XXZZ,XXYZ,XXYY,XXXZ,XXXY,XXXX
	!
	!G 0=ZZZZ+3/8*(XXXX+YYYY)-3*√(3/35)*(XXZZ+YYZZ-1/4*XXYY)
	 matg(1,1)=1D0
	 matg(3,1)=-3D0*sqrt(3D0/35D0)
	 matg(5,1)=3D0/8D0
	 matg(10,1)=-3D0*sqrt(3D0/35D0)
	 matg(12,1)=3D0/4D0*sqrt(3D0/35D0)
	 matg(15,1)=3D0/8D0
	 !G+1=2*√(5/14)*XZZZ-3/2*√(5/14)*XXXZ-3/2/√14*XYYZ
	 matg(6,2)=2D0*sqrt(5D0/14D0)
	 matg(8,2)=-1.5D0/sqrt(14D0)
	 matg(13,2)=-1.5D0*sqrt(5D0/14D0)
	 !G-1=2*√(5/14)*YZZZ-3/2*√(5/14)*YYYZ-3/2/√14*XXYZ
	 matg(2,3)=2D0*sqrt(5D0/14D0)
	 matg(4,3)=-1.5D0*sqrt(5D0/14D0)
	 matg(11,3)=-1.5D0/sqrt(14D0)
	 !G+2=3*√(3/28)*(XXZZ-YYZZ)-√5/4*(XXXX-YYYY)
	 matg(3,4)=-3D0*sqrt(3D0/28D0)
	 matg(5,4)=sqrt(5D0)/4D0
	 matg(10,4)=3D0*sqrt(3D0/28D0)
	 matg(15,4)=-sqrt(5D0)/4D0
	 !G-2=3/√7*XYZZ-√(5/28)*(XXXY+XYYY)
	 matg(7,5)=3D0/sqrt(7D0)
	 matg(9,5)=-sqrt(5D0/28D0)
	 matg(14,5)=-sqrt(5D0/28D0)
	 !G+3=√(5/8)*XXXZ-3/√8*XYYZ
	 matg(8,6)=-3D0/sqrt(8D0)
	 matg(13,6)=sqrt(5D0/8D0)
	 !G-3=-√(5/8)*YYYZ+3/√8*XXYZ
	 matg(4,7)=-sqrt(5D0/8D0)
	 matg(11,7)=3D0/sqrt(8D0)
	 !G+4=√35/8*(XXXX+YYYY)-3/4*√3*XXYY
	 matg(5,8)=sqrt(35D0)/8D0
	 matg(12,8)=-3D0/4D0*sqrt(3D0)
	 matg(15,8)=sqrt(35D0)/8D0
	 !G-4=√5/2*(XXXY-XYYY)
	 matg(9,9)=-sqrt(5D0)/2D0
	 matg(14,9)=sqrt(5D0)/2D0
else if (iprog==2) then !For .molden
	! From 9G: G 0,G+1,G-1,G+2,G-2,G+3,G-3,G+4,G-4
	! To 15G:   1    2    3    4    5    6    7    8
	!         xxxx,yyyy,zzzz,xxxy,xxxz,yyyx,yyyz,zzzx
	!           9   10   11   12   13   14   15
	!         zzzy,xxyy,xxzz,yyzz,xxyz,yyxz,zzxy
	!
	!G 0=ZZZZ+3/8*(XXXX+YYYY)-3*√(3/35)*(XXZZ+YYZZ-1/4*XXYY)
	matg(3,1)=1D0
	matg(1,1)=3D0/8D0
	matg(2,1)=3D0/8D0
	matg(11,1)=-3D0*sqrt(3D0/35D0)
	matg(12,1)=-3D0*sqrt(3D0/35D0)
	matg(10,1)=3D0/4D0*sqrt(3D0/35D0)
	!G+1=2*√(5/14)*XZZZ-3/2*√(5/14)*XXXZ-3/2/√14*XYYZ
	matg(8,2)=2D0*sqrt(5D0/14D0)
	matg(5,2)=-1.5D0*sqrt(5D0/14D0)
	matg(14,2)=-1.5D0/sqrt(14D0)
	!G-1=2*√(5/14)*YZZZ-3/2*√(5/14)*YYYZ-3/2/√14*XXYZ
	matg(9,3)=2D0*sqrt(5D0/14D0)
	matg(7,3)=-1.5D0*sqrt(5D0/14D0)
	matg(13,3)=-1.5D0/sqrt(14D0)
	!G+2=3*√(3/28)*(XXZZ-YYZZ)-√5/4*(XXXX-YYYY)
	matg(11,4)=3D0*sqrt(3D0/28D0)
	matg(12,4)=-3D0*sqrt(3D0/28D0)
	matg(1,4)=-sqrt(5D0)/4D0
	matg(2,4)=sqrt(5D0)/4D0
	!G-2=3/√7*XYZZ-√(5/28)*(XXXY+XYYY)
	matg(15,5)=3D0/sqrt(7D0)
	matg(4,5)=-sqrt(5D0/28D0)
	matg(6,5)=-sqrt(5D0/28D0)
	!G+3=√(5/8)*XXXZ-3/√8*XYYZ
	matg(5,6)=sqrt(5D0/8D0)
	matg(14,6)=-3D0/sqrt(8D0)
	!G-3=-√(5/8)*YYYZ+3/√8*XXYZ
	matg(7,7)=-sqrt(5D0/8D0)
	matg(13,7)=3D0/sqrt(8D0)
	!G+4=√35/8*(XXXX+YYYY)-3/4*√3*XXYY
	matg(1,8)=sqrt(35D0)/8D0
	matg(2,8)=sqrt(35D0)/8D0
	matg(10,8)=-3D0/4D0*sqrt(3D0)
	!G-4=√5/2*(XXXY-XYYY)
	matg(4,9)=sqrt(5D0)/2D0
	matg(6,9)=-sqrt(5D0)/2D0
end if

! From 11H: H 0,H+1,H-1,H+2,H-2,H+3,H-3,H+4,H-4,H+5,H-5
! To 21H:   1     2     3     4     5     6     7     8     9    10
!         ZZZZZ YZZZZ YYZZZ YYYZZ YYYYZ YYYYY XZZZZ XYZZZ XYYZZ XYYYZ 
!          11    12    13    14    15    16    17    18    19    20    21
!         XYYYY XXZZZ XXYZZ XXYYZ XXYYY XXXZZ XXXYZ XXXYY XXXXZ XXXXY XXXXX
!
!H 0=ZZZZZ-5/√21*(XXZZZ+YYZZZ)+5/8*(XXXXZ+YYYYZ)+√(15/7)/4*XXYYZ
math(1,1)=1D0
math(12,1)=-5D0/sqrt(21D0)
math(3,1)=-5D0/sqrt(21D0)
math(19,1)=5D0/8D0
math(5,1)=5D0/8D0
math(14,1)=sqrt(15D0/7D0)/4D0
!H+1=√(5/3)*XZZZZ-3*√(5/28)*XXXZZ-3/√28*XYYZZ+√15/8*XXXXX+√(5/3)/8*XYYYY+√(5/7)/4*XXXYY
math(7,2)=sqrt(5D0/3D0)
math(16,2)=-3D0*sqrt(5D0/28D0)
math(9,2)=-3D0/sqrt(28D0)
math(21,2)=sqrt(15D0)/8D0
math(11,2)=sqrt(5D0/3D0)/8D0
math(18,2)=sqrt(5D0/7D0)/4D0
!H-1=√(5/3)*YZZZZ-3*√(5/28)*YYYZZ-3/√28*XXYZZ+√15/8*YYYYY+√(5/3)/8*XXXXY+√(5/7)/4*XXYYY
math(2,3)=sqrt(5D0/3D0)
math(4,3)=-3D0*sqrt(5D0/28D0)
math(13,3)=-3D0/sqrt(28D0)
math(6,3)=sqrt(15D0)/8D0
math(20,3)=sqrt(5D0/3D0)/8D0
math(15,3)=sqrt(5D0/7D0)/4D0
!H+2=√5/2*(XXZZZ-YYZZZ)-√(35/3)/4*(XXXXZ-YYYYZ)
math(12,4)=sqrt(5D0)/2D0
math(3,4)=-sqrt(5D0)/2D0
math(19,4)=-sqrt(35D0/3D0)/4D0
math(5,4)=sqrt(35D0/3D0)/4D0
!H-2=√(5/3)*XYZZZ-√(5/12)*(XXXYZ+XYYYZ)
math(8,5)=sqrt(5D0/3D0)
math(17,5)=-sqrt(5D0/12D0)
math(10,5)=-sqrt(5D0/12D0)
!H+3=√(5/6)*XXXZZ-√(3/2)*XYYZZ-√(35/2)/8*(XXXXX-XYYYY)+√(5/6)/4*XXXYY
math(16,6)=sqrt(5D0/6D0)
math(9,6)=-sqrt(1.5D0)
math(21,6)=-sqrt(17.5D0)/8D0
math(11,6)=sqrt(17.5D0)/8D0
math(18,6)=sqrt(5D0/6D0)/4D0
!H-3=-√(5/6)*YYYZZ+√(3/2)*XXYZZ-√(35/2)/8*(XXXXY-YYYYY)-√(5/6)/4*XXYYY
math(4,7)=-sqrt(5D0/6D0)
math(13,7)=sqrt(1.5D0)
math(20,7)=-sqrt(17.5D0)/8D0
math(6,7)=sqrt(17.5D0)/8D0
math(15,7)=-sqrt(5D0/6D0)/4D0
!H+4=√35/8*(XXXXZ+YYYYZ)-3/4*√3*XXYYZ
math(19,8)=sqrt(35D0)/8D0
math(5,8)=sqrt(35D0)/8D0
math(14,8)=-0.75D0*sqrt(3D0)
!H-4=√5/2*(XXXYZ-XYYYZ)
math(17,9)=sqrt(5D0)/2D0
math(10,9)=-sqrt(5D0)/2D0
!H+5=3/8*√(7/2)*XXXXX+5/8*√(7/2)*XYYYY-5/4*√(3/2)*XXXYY
math(21,10)=3D0/8D0*sqrt(3.5D0)
math(11,10)=5D0/8D0*sqrt(3.5D0)
math(18,10)=-1.25D0*sqrt(1.5D0)
!H-5=3/8*√(7/2)*YYYYY+5/8*√(7/2)*XXXXY-5/4*√(3/2)*XXYYY
math(6,11)=3D0/8D0*sqrt(3.5D0)
math(20,11)=5D0/8D0*sqrt(3.5D0)
math(15,11)=-1.25D0*sqrt(1.5D0)
end subroutine




!!---------- Load Fock or Kohn-Sham matrix from NBO .47 file, ORCA output file or plain text file
!User can also require to directly generate Fock matirx based on orbital energies
!istatus=0 means successfully loaded. =1 means failed
subroutine loadFockfile(istatus)
use defvar
use util
character c200tmp*200
integer istatus

do while(.true.)
	write(*,"(/,a)") " Input the file recording Fock/KS matrix in original basis functions in lower triangular form, e.g. C:\Fock.txt"
	write(*,"(a)") " Note: If file suffix is .47 or .mwfn, then Fock/KS matrix will be directly loaded from it"
	write(*,"(a)") "       If this is an ORCA output file using ""%output Print[P_Iter_F] 1 end"", Fock/KS matrix printed at last iteration will be loaded"
    !Not quite meaningful, because when need to use Fock matrix, current orbitals are not original MOs
	!write(*,*) "To directly generate Fock/KS matrix based on orbital energies, input ""gen"""
	read(*,"(a)") c200tmp
    if (c200tmp=="gen") then
		call MOene2Fmat(istatus)
        if (istatus==0) return !Successfully generated
    else
		inquire(file=c200tmp,exist=alive)
		if (.not.alive) then
			write(*,*) "Error: Unable to find this file!"
			cycle
		end if
		exit
    end if
end do

open(10,file=c200tmp,status="old")
call outputprog(10,iprog)
if (allocated(FmatA)) deallocate(FmatA)
allocate(FmatA(nbasis,nbasis))
if (iprog==2) then !ORCA output file
    write(*,*) "This file is recognized as an ORCA output file"
    call loclabelfinal(10,"Fock matrix for operator 0",nfound)
    if (nfound==0) then
		write(*,*) "Error: Unable to locate ""Fock matrix for operator 0"" in this file!"
		close(10)
		istatus=1
		return
    end if
    if (wfntype==1.or.wfntype==4) then !Unrestricted
		write(*,*) "Loading alpha Fock/KS matrix (Fock matrix for operator 0)..."
    else
		write(*,*) "Loading Fock/KS matrix (Fock matrix for operator 0)..."
    end if
    call readmatgau(10,FmatA,0,"?",10,6)
    call ORCAmatconv(nbasis,FmatA)
else !.47 or .mwfn .or. plain text file
    if (index(c200tmp,".47")/=0) then
	    if (wfntype==1.or.wfntype==4) then !Unrestricted
			write(*,*) "Trying to load alpha Fock/KS matrix from .47 file..."
	    else
			write(*,*) "Trying to load Fock/KS matrix from .47 file..."
        end if
	    call loclabel(10,"$FOCK",ifound)
	    if (ifound==0) then
		    write(*,*) "Error: Unable to find $FOCK field in this file!"
		    close(10)
		    istatus=1
		    return
	    end if
	    read(10,*)
    else if (index(c200tmp,".mwfn")/=0) then
	    if (wfntype==1.or.wfntype==4) then !Unrestricted
			write(*,*) "Trying to load alpha Fock/KS matrix from .mwfn file..."
	    else
			write(*,*) "Trying to load Fock/KS matrix from .mwfn file..."
        end if
	    call loclabel(10,"1-e Hamiltonian matrix",ifound) !For R/RO case, locate to $1-e Hamiltonian matrix, for U case, locate to $Alpha 1-e Hamiltonian matrix
	    if (ifound==0) then
		    write(*,*) "Error: Unable to find 1-e Hamiltonian matrix in this file!"
		    close(10)
		    istatus=1
		    return
	    end if
	    read(10,*)
    else !Plain text file
	    rewind(10)
    end if
    read(10,*) ((FmatA(i,j),j=1,i),i=1,nbasis) !Load total or alpha Fock matrix (lower triangular part)
    do i=1,nbasis !Fill upper triangular part
	    do j=i+1,nbasis
		    FmatA(i,j)=FmatA(j,i)
	    end do
    end do
end if
!Checking detail of loaded Fock matrix
!call showmatgau(FmatA,form="f14.6")
!ibas=0
!do ish=1,nshell
!    ishtype=shtype(ish)
!    nshbas=shtype2nbas(ishtype)
!    if (shtype(ish)==-3) then
!        call showmatgau(FmatA(:,ibas+1:ibas+nshbas),form="f14.6")
!        exit
!    end if
!    ibas=ibas+nshbas
!end do
if (wfntype==1.or.wfntype==4) then !Also load beta part
	if (allocated(FmatB)) deallocate(FmatB)
	allocate(FmatB(nbasis,nbasis))
    if (iprog==2) then !ORCA
        write(*,*) "Loading beta Fock/KS matrix (Fock matrix for operator 1)..."
        call readmatgau(10,FmatB,0,"?",10,6)
        call ORCAmatconv(nbasis,FmatB)
    else
	    write(*,*) "Trying to load beta Fock/KS matrix from the file..."
        if (index(c200tmp,".mwfn")/=0) then
            call loclabel(10,"Beta 1-e Hamiltonian matrix",ifound)
            read(10,*)
        end if
	    read(10,*,iostat=ierror) ((FmatB(i,j),j=1,i),i=1,nbasis) !Load beta Fock matrix
        if (ierror==0) then
			do i=1,nbasis !Fill upper triangular part
				do j=i+1,nbasis
					FmatB(i,j)=FmatB(j,i)
				end do
			end do
        else
			write(*,"(a)") " Warning: Unable to load Fock/KS matrix of beta spin! Now the beta Fock/KS matrix has been set to alpha Fock/KS instead"
            write(*,*) "Press ENTER button to continue"
            read(*,*)
            FmatB=FmatA
        end if
    end if
end if
close(10)
write(*,*) "Fock/KS matrix has been loaded successfully!"
istatus=0
end subroutine



!!-------- Randomly generate name of Sobereva's lover
subroutine mylover(outname)
integer,parameter :: nlovers=59
character(len=80) lovername(nlovers),outname
CALL RANDOM_SEED()
CALL RANDOM_NUMBER(tmp)
lovername(1)="K-ON\Mio_Akiyama"
lovername(2)="K-ON\Azusa_Nakano"
lovername(3)="EVA\Rei_Ayanami"
lovername(4)="Ore_no_Imoto\Black_Cat"
lovername(5)="Touhou_project\Ran_Yakumo"
lovername(6)="Haiyore!Nyaruko-san\Nyaruko"
lovername(7)="Bodacious_Space_Pirates\Kurihara_Chiaki"
lovername(8)="Otoboku\Mariya_Mikado"
lovername(9)="Amagami\Miya_Tachibana"
lovername(10)="Shakugan_no_Shana\Shana"
lovername(11)="Tiger_Mask_W\Miss_X"
! lovername(11)="Yuru_Yuri\Akari_Akaza"
lovername(12)="Natsuiro_Kiseki\Yuka_Hanaki"
lovername(13)="Love_Live!\Nico_Yazawa"
lovername(14)="Love_Live!\Nozomi_Tojo"
lovername(15)="Love_Live!\Nishikino_Maki"
lovername(16)="Last_Exile\Dio_Eraclea"
lovername(17)="NHK_ni_Youkoso!\Misaki_Nakahara"
lovername(18)="Rio_Rainbow_Gate\Rio_Rollins"
lovername(19)="Blood-C\Saya_Kisaragi"
lovername(20)="Mahou_Shoujo_Madoka-Magica\Homura_Akemi"
lovername(21)="Saki\Hisa_Takei"
lovername(22)="Strawberry_Panic\Chikaru_Minamoto"
lovername(23)="Najica\Najica_Hiiragi"
lovername(24)="Blue_Drop\Hagino_Senkouji"
lovername(25)="Fate_Zero\Saber"
lovername(26)="Baka_to_Test_to_Shoukanjuu\Hideyoshi_Kinoshita"
lovername(27)="Watamote\Tomoko_Kuroki"
lovername(28)="Genshiken_Nidaime\Kenjirou_Hato"
lovername(29)="Love_is_Like_After_the_Rain\Akira_Tachibana"
lovername(30)="Kan_Colle\Shimakaze"
lovername(31)="Kan_Colle\Kongou"
lovername(32)="Gokukoku\Kazumi_Schlierenzauer"
lovername(33)="Vocaloid\Miku_Hatsune"
lovername(34)="Tokimeki_Memorial\Yuina_Himoo"
lovername(35)="MADLAX\MADLAX"
lovername(36)="Gun_Gale_Online\Kirito"
lovername(37)="Denkigai_No_Honyasan\Sennsei"
lovername(38)="Wake_Up,Girls!\Miyu_Okamoto"
lovername(39)="Plastic_Memories\Aira"
lovername(40)="Real_world\Kanan" !"Real_world\sell-moe-kun"
lovername(41)="Sakurako-san_no_Ashimoto_ni_wa_Shitai_ga_Umatteiru\Sakurako"
lovername(42)="Hibike!_Euphonium\Reina_Kousaka"
lovername(43)="Planetarian\Yumemi_Hoshino"
lovername(44)="Lovelive_Sunshine!!\Yoshiko_Tsushima"
lovername(45)="Lovelive_Sunshine!!\Kurosawa Dia"
lovername(46)="Lovelive_Sunshine!!\Sakurauchi_Riko"
lovername(47)="Violet_Evergarden\Violet_Evergarden"
lovername(48)="Otobuko\Mizuho_Miyanokouji"
lovername(49)="iDOLM@STER\Makoto_Kikuchi"
lovername(50)="Fate\Rin_Tohsaka"
lovername(51)="Magical_Girl_Spec-Ops_Asuka\Asuka_Otori"
lovername(52)="Granblue_Fantasy\Katalina"
lovername(53)="Gochuumon_wa_Usagi_Desu_ka\Rize_Tedeza"
lovername(54)="Date_A_Live\Kurumi_Tokisaki"
lovername(55)="Adachi_to_shimamura\Sakura_Adachi"
lovername(56)="Don't_Toy_with_Me,_Miss_Nagatoro\Hayase_Nagatoro"
lovername(57)="Super_Cub\Reiko"
lovername(58)="LoveLive!_Superstar!!\Sumire_Heanna"
lovername(59)="Jahy-sama_wa_Kujikenai!\Jahy"

!Dear Kanan,
!
!You are the only one I deeply love forever in the real world,
!although you can't be with me, and I am even unable to know your name and touch your finger.
!I believe I will never love anyone else in the rest of my life.
!
!I love your brilliant dance, your kawaii smile, your lovely double ponytail, and especially, your extremely pure and beautiful heart.
!
!                     ----- 2015-May-19
outname=lovername(ceiling(tmp*nlovers))
end subroutine





!!----------- Convert current CObasa / CObasb to CO
!ispin=1: Only alpha, =2: Only beta, =3: Both alpha and beta
subroutine CObas2CO(ispin)
use defvar
implicit real*8 (a-h,o-z)
integer ispin
real*8 conv5d6d(6,5),conv7f10f(10,7),conv9g15g(15,9),conv11h21h(21,11)
real*8,allocatable :: CObasa_cart(:,:),CObasb_cart(:,:)

call gensphcartab(1,conv5d6d,conv7f10f,conv9g15g,conv11h21h)

nbasis_cart=sum(shtype2nbas(abs(shtype(:))))
if (ispin==1.or.ispin==3) allocate(CObasa_cart(nbasis_cart,nbasis))
if (ispin==2.or.ispin==3) allocate(CObasb_cart(nbasis_cart,nbasis))
CObasa_cart=0

!Map spherical coefficients to Cartesian coefficients
ipos5D=1
ipos6D=1
do ish=1,nshell
    ishtype=shtype(ish)
    numshbas5D=shtype2nbas(ishtype)
    numshbas6D=shtype2nbas(abs(ishtype))
    if (ispin==1.or.ispin==3) then !Alpha part
        if (ishtype>=0) then !S,P or Cartesian type, in this case numshbas5D=numshbas6D 
            CObasa_cart(ipos6D:ipos6D+numshbas6D-1,:)=CObasa(ipos5D:ipos5D+numshbas5D-1,:)
        else
	        if (ishtype==-2) then
		        CObasa_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv5d6d,CObasa(ipos5D:ipos5D+numshbas5D-1,:))
	        else if (ishtype==-3) then
		        CObasa_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv7f10f,CObasa(ipos5D:ipos5D+numshbas5D-1,:))
	        else if (ishtype==-4) then
		        CObasa_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv9g15g,CObasa(ipos5D:ipos5D+numshbas5D-1,:))
	        else if (ishtype==-5) then
		        CObasa_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv11h21h,CObasa(ipos5D:ipos5D+numshbas5D-1,:))
	        end if
        end if
    end if
    if (ispin==2.or.ispin==3) then !Beta part
        if (ishtype>=0) then !S,P or Cartesian type, in this case numshbas5D=numshbas6D 
            CObasb_cart(ipos6D:ipos6D+numshbas6D-1,:)=CObasb(ipos5D:ipos5D+numshbas5D-1,:)
        else
	        if (ishtype==-2) then
		        CObasb_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv5d6d,CObasb(ipos5D:ipos5D+numshbas5D-1,:))
	        else if (ishtype==-3) then
		        CObasb_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv7f10f,CObasb(ipos5D:ipos5D+numshbas5D-1,:))
	        else if (ishtype==-4) then
		        CObasb_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv9g15g,CObasb(ipos5D:ipos5D+numshbas5D-1,:))
	        else if (ishtype==-5) then
		        CObasb_cart(ipos6D:ipos6D+numshbas6D-1,:)=matmul(conv11h21h,CObasb(ipos5D:ipos5D+numshbas5D-1,:))
	        end if
        end if
    end if
	ipos5D=ipos5D+numshbas5D
	ipos6D=ipos6D+numshbas6D
end do

do imo=1,nbasis
    do ibas=1,nbasis_cart
        if (ispin==1.or.ispin==3) then
            do iGTF=primstart(ibas),primend(ibas)
                CO(imo,iGTF)=CObasa_cart(ibas,imo)*primconnorm(iGTF)
            end do
        end if
        if (ispin==2.or.ispin==3) then
            do iGTF=primstart(ibas),primend(ibas)
                CO(imo+nbasis,iGTF)=CObasb_cart(ibas,imo)*primconnorm(iGTF)
            end do
        end if
    end do
end do

end subroutine




!!------- Add a Bq atom to specific position
subroutine addBq(xpos,ypos,zpos)
use defvar
real*8 xpos,ypos,zpos
if (allocated(a_tmp)) deallocate(a_tmp)
allocate(a_tmp(ncenter))
a_tmp=a
ncenter=ncenter+1
deallocate(a)
allocate(a(ncenter))
a(1:ncenter-1)=a_tmp
a(ncenter)%index=0
a(ncenter)%charge=0
a(ncenter)%name="Bq"
a(ncenter)%x=xpos
a(ncenter)%y=ypos
a(ncenter)%z=zpos
deallocate(a_tmp)
end subroutine


!!------- Invoke Gaussian to run a .gjf
!If returned istate=1, means normally termination, =0 means other case or failed
subroutine runGaussian(gjfname,istate)
use defvar
use util
character(len=*) gjfname
character command*400,outname*200
outname=gjfname(:len(gjfname)-3)//"out"
command='"'//trim(gaupath)//'" "'//gjfname//'" "'//trim(outname)//'"'
call runcommand(command)
open(100,file=outname,status="old")
call loclabel(100,"Normal termination",istate)
close(100)
end subroutine


!!------- Invoke ORCA to run a .inp
!If returned istate=1, means normally termination, =0 means other case or failed
subroutine runORCA(inpname,istate)
use defvar
use util
character(len=*) inpname
character command*400,outname*200
outname=inpname(:len(inpname)-3)//"out"
command='"'//trim(orcapath)//'" "'//inpname//'" > "'//trim(outname)//'"'
call runcommand(command)
open(100,file=outname,status="old")
call loclabel(100,"****ORCA TERMINATED NORMALLY****",istate)
close(100)
end subroutine



!!---------- Generate connectivity matrix
!infomode =1: Output prompts  =0: Silent
!iallowPBC =1: Allow considering PBC when this system is PBC  =0: Not allow (can avoid incompatibility problem)
subroutine genconnmat(infomode,iallowPBC)
use defvar
use util
implicit real*8 (a-h,o-z)
integer infomode,iallowPBC

!call walltime(iwalltime1)
if (allocated(a)) then
    if (allocated(connmat)) deallocate(connmat)
    allocate(connmat(ncenter,ncenter))
    if (infomode==1) then
		write(*,*) "Generating bonding relationship..."
		write(*,"(a,f5.3,a)") " Note: If distance between two atoms is smaller than sum of their &
		covalent radii multiplied by ",bondcrit,", then they are regarded as bonded"
    end if
    connmat=0
	!$OMP PARALLEL DO SHARED(connmat) PRIVATE(iatm,jatm,covri,covrj,tmpdist) schedule(dynamic) NUM_THREADS(nthreads)
    do iatm=1,ncenter
		covri=covr(a(iatm)%index)
        do jatm=iatm+1,ncenter
			if (iallowPBC==1.and.ifPBC>0) then
				call nearest_atmdistxyz(iatm,jatm,tmpdist,atmx,atmy,atmz)
            else	
				tmpdist=atomdist(iatm,jatm)
            end if
			covrj=covr(a(jatm)%index)
            if ( tmpdist < bondcrit*(covri+covrj) ) connmat(iatm,jatm)=1
            connmat(jatm,iatm)=connmat(iatm,jatm)
        end do
    end do
	!$OMP END PARALLEL DO
else
    if (infomode==1) write(*,"(a)") " Unable to generate bonding relationship because there is no atom information!"
end if

call walltime(iwalltime2)
!write(*,"(' Generating connectivity matrix took up wall clock time',i10,' s')") iwalltime2-iwalltime1

!Check connectivity
!do i=1,ncenter
!    do j=i+1,ncenter
!        if (connmat(i,j)==1) write(*,*) i,j,connmat(i,j)
!    end do
!end do
end subroutine



!!--------- Generate neighbouring list of each atom
!itype=1: Based on existing connectivity matrix "connmat"
!itype=2: Directly based on interatomic distance
subroutine genneighlist(itype,iallowPBC,neigh,nneigh)
use defvar
use util
implicit real*8 (a-h,o-z)
integer nneigh(ncenter) !Number of neighbours (the number of atoms connected to this atom)
integer neigh(maxneigh,ncenter) !neigh(1:nneigh(i),i) is list of neighbouring atom indices of atom i
integer itype,iallowPBC

nneigh=0
if (itype==1) then
	if (.not.allocated(connmat)) call genconnmat(1,iallowPBC) !Generate connectivity matrix
	!$OMP PARALLEL DO SHARED(nneigh,neigh) PRIVATE(iatm,jatm) schedule(dynamic) NUM_THREADS(nthreads)
	do iatm=1,ncenter
		do jatm=1,ncenter
			if (connmat(jatm,iatm)>0) then
				nneigh(iatm)=nneigh(iatm)+1
				neigh(nneigh(iatm),iatm)=jatm
			end if
		end do
	end do
	!$OMP END PARALLEL DO
else if (itype==2) then
	write(*,*) "Generating bonding relationship..."
	write(*,"(a,f5.3,a)") " Note: If distance between two atoms is smaller than sum of their &
	covalent radii multiplied by ",bondcrit,", then they are regarded as bonded"
    !$OMP PARALLEL DO SHARED(nneigh,neigh) PRIVATE(iatm,jatm,covri,covrj,tmpdist) schedule(dynamic) NUM_THREADS(nthreads)
    do iatm=1,ncenter
		covri=covr(a(iatm)%index)
        do jatm=1,ncenter
			if (iallowPBC==1.and.ifPBC>0) then
				call nearest_atmdistxyz(iatm,jatm,tmpdist,atmx,atmy,atmz)
            else	
				tmpdist=atomdist(iatm,jatm)
            end if
			covrj=covr(a(jatm)%index)
            if ( tmpdist < bondcrit*(covri+covrj) ) then
				nneigh(iatm)=nneigh(iatm)+1
				neigh(nneigh(iatm),iatm)=jatm
			end if
        end do
    end do
	!$OMP END PARALLEL DO
end if

!write(*,*) "Neighbouring list:"
!do iatm=1,ncenter
!	write(*,"(i6,a,':')",advance="no") iatm,a(iatm)%name
!	do idx=1,nneigh(iatm)
!		write(*,"(i6)",advance="no") neigh(idx,iatm)
!    end do
!    write(*,*)
!end do
end subroutine



!!--------- Generate fragment index of all atoms according to connectivity
!Each atom has initial fragment index of itself, then each atom is compare with adjacent ones, if an adjacent atom has fragment &
!index smaller than it, the its fragment will be updated to the adjacent one. Iteration performs until no fragment updation is occurred
subroutine genconnfrag(atmfrg)
use defvar
use util
integer atmfrg(ncenter),tmpidx(ncenter)
integer neigh(maxneigh,ncenter),nneigh(ncenter)
character c2000tmp*2000

write(*,*) "Generate neighbouring list..."
call genneighlist(2,1,neigh,nneigh) !Generate according to interatomic distance

!write(*,*) "Generate atmfrg..."
forall (i=1:ncenter) atmfrg(i)=i
do while(.true.)
    inew=0
    do iatm=1,ncenter
        do jdx=1,nneigh(iatm)
			jatm=neigh(jdx,iatm)
            if (atmfrg(jatm)<atmfrg(iatm)) then
                atmfrg(iatm)=atmfrg(jatm)
                inew=inew+1
                exit
            end if
        end do
    end do
    if (inew==0) exit
end do

!Make fragment indices contiguous
nfrg=0
do iatm=1,ncenter
    natmfrg=count(atmfrg==iatm)
    if (natmfrg/=0) then
        nfrg=nfrg+1
        where(atmfrg==iatm) atmfrg=nfrg
    end if
end do

!Show fragment information
!do ifrg=1,nfrg
!    write(*,"(' Fragment',i7,', number of atoms',i7)") ifrg,count(atmfrg==ifrg)
!    itmp=0
!    do iatm=1,ncenter
!        if (atmfrg(iatm)==ifrg) then
!            itmp=itmp+1
!            tmpidx(itmp)=iatm
!        end if
!    end do
!    call arr2str_2(tmpidx(1:itmp),c2000tmp)
!    write(*,"(a)") trim(c2000tmp)
!end do
end subroutine



!!-------- Input index of an atom, then the indices of all atoms in the fragment will be returned
!"iatm" is the selected atom, "iffrag" has length of ncenter, if an atom is in the fragment, the value is 1, else 0
subroutine getfragatoms(iselatm,iffrag)
use defvar
use util
implicit real*8 (a-h,o-z)
integer iselatm,iffrag(ncenter)
iffrag=0
iffrag(iselatm)=1
if (.not.allocated(connmat)) call genconnmat(0,0) !Generate connectivity matrix
do while(.true.)
    inew=0
    do iatm=1,ncenter !Cycle all atoms, if it is not in fragment, and it is linked to an atom already in fragment, it will be added to fragment
        if (iffrag(iatm)==1) cycle !The atom is already in fragment
        if (ishowhydrogen==0.and.a(iatm)%index==1) cycle !If this is hydrogen while we request do not showing hydrogens, skip it
        do jatm=1,ncenter !Cycle neighbouring atoms
            if (jatm==iatm) cycle
            if (connmat(iatm,jatm)>0.and.iffrag(jatm)==1) then
                iffrag(iatm)=1
                inew=inew+1
                exit
            end if
        end do
    end do
    if (inew==0) exit
end do
end subroutine



!!--------- Align atoms in "a" (mol1) to "a_org" (mol2), the atoms must have the same order and same number
!See http://nghiaho.com/?page_id=671 "FINDING OPTIMAL ROTATION AND TRANSLATION BETWEEN CORRESPONDING 3D POINTS"
!Test purpose:
!call readinfile("C:\Users\Sobereva\Desktop\md2.xyz",0)
!call geomalign
!call outxyz("new.xyz",10)
subroutine geomalign
use defvar
use util
implicit real*8 (a-h,o-z)
real*8 xyz1(3,ncenter),xyz2(3,ncenter),eleidx(ncenter),com1(3),com2(3)  !1/2/3=x,y,z
real*8 singval(3),Umat(3,3),Vmat(3,3),Hmat(3,3),Rmat(3,3)

do i=1,ncenter
    xyz1(1,i)=a_org(i)%x
    xyz1(2,i)=a_org(i)%y
    xyz1(3,i)=a_org(i)%z
    xyz2(1,i)=a(i)%x
    xyz2(2,i)=a(i)%y
    xyz2(3,i)=a(i)%z
    eleidx(i)=a_org(i)%index
end do
totmass=sum(atmwei(eleidx(:)))
do i=1,3
    com1(i)=sum( atmwei(eleidx(:))*xyz1(i,:) )/totmass
    com2(i)=sum( atmwei(eleidx(:))*xyz2(i,:) )/totmass
end do

Hmat=0
do iatm=1,ncenter
    xyz1(:,iatm)=xyz1(:,iatm)-com1(:)
    xyz2(:,iatm)=xyz2(:,iatm)-com2(:)
    Hmat=Hmat+matmul(xyz2(:,iatm:iatm),transpose(xyz1(:,iatm:iatm)))
end do

call SVDmat(1,Hmat,Umat,Vmat,singval,info)
Rmat=matmul(Vmat,transpose(Umat))
!write(*,*) detmat(Rmat)
!call showmatgau(Rmat,"Rotation matrix")
!if (detmat(Rmat)<0) Rmat(:,3)=-Rmat(:,3) !As mentioned in the webpage, this line should be added, however it incorrectly make molecule mirror inverted!

!Rotate mol2
do i=1,ncenter
    xyz2(:,i)=matmul(Rmat,xyz2(:,i))
end do

!Move center of mass of mol2 to mol1
do i=1,ncenter
    xyz2(:,i)=xyz2(:,i)+com1(:)
end do

!RMSD
RMSD=0
do i=1,ncenter
    RMSD=RMSD+sum(((xyz1(:,i)+com1(:))-xyz2(:,i))**2)
end do
RMSD=dsqrt(RMSD/ncenter)
write(*,"(' RMSD:',f12.6,' Bohr')") RMSD

do i=1,ncenter
    a(i)%x=xyz2(1,i)
    a(i)%y=xyz2(2,i)
    a(i)%z=xyz2(3,i)
end do
end subroutine



!!------- Get point group and list of symmetry-equivalence atoms by invoking SYVA routines
!The tolerance affects if symmetry-equivalence atoms and point group could be successfully recognized, the value
!should not be too large, otherwise the recognition may be completely failed!
!  Input variables:
!natoms: The number of inputted atoms
!nat(1:natoms): Element index of all inputted atoms
!coord(1:3,1:natoms): x,y,z of all inputted atoms, must be in Angstrom, otherwise point group determination may be incorrect!
!delta: Representing the distortion of the geometry. 0.01 is aproximately identical to "default" in gview, 0.1 corresponds to loose
!  Returned variables:
!pglabel: Point group label
!nclass: The number of equivalent classes, should have size of "natoms"
!classnatm(i): The number of atoms in equivalent class i, should have size of "natoms"
!classidx(:,i): The atom indices in equivalent class i, both dimension should have size of "natoms"
subroutine PG_eqvatm(natoms,nat,coord,delta,pglabel,nclass,classnatm,classidx)
implicit real*8 (a-h,o-z)
common /data/ wt(90),symb(90)
common /chartab/ nir(2,55),chtab(14,322),nsymop(14,4,55),nrotharm(3,322),pgsymb(57),irsymb(322)
common /subgroups/ nsgb(2,57),nsgr(406)
character symb*2,pglabel*3,pgsymb*3,irsymb*4
integer natoms,nclass
integer classnatm(natoms),classidx(natoms,natoms)
integer,parameter :: nmax=200 !nmat: maximum number of symmetic operation
real*8 coord(3,natoms),delta,pc(3),symn(3,nmax)
integer nat(natoms),nper(natoms,250),nscl(natoms,natoms),nccl(natoms),nsym(nmax,5)
ncr=0
nsr=0
nsg=0
nout=0 !Suppress almost all output of SYVA routines

if (natoms==1) then
    nclass=1
    classnatm(1)=1
    classidx(1,1)=1
    pglabel="C1"
    return
end if

!Check atoms passed-in
!write(*,'(1x,a,10x,a1,16x,a1,15x,a1,15x,a1)') 'number','x','y','z','w'
!do i=1,natoms
!   write(*,'(1x,i3,2x,4f16.6)') nat(i),(coord(j,i),j=1,3),wt(nat(i))
!end do

!Calculation of the COM (centre of mass) of the molecule
call syva_cmass(natoms,nat,wt,coord,wmol,cmx,cmy,cmz)
pc(1)=cmx;pc(2)=cmy;pc(3)=cmz

!Shift the origin of the Cartesian system to COM
call syva_cshift(natoms,coord,pc)

!Find symmetry operations
call sym_elements(natoms,nat,coord,symb,delta,ng,ni,nsg,ncr,nsr, np,symn,nsym,nout,nprm,nper,nseq,nccl,nscl)

!Detemines the equivalence classes defined by the symmetry operations
call symclass(natoms,nprm,nper,nseq,nccl,nscl,nat,symb,nout)

!write(*,*) "Symmetry-equivalence classes of atoms: ",nseq
!do i=1,nseq
!   write(*,'(/5x,a,i3,a7,a2,a1)') '#',i,' (atom ',symb(nat(nscl(1,i))),')'
!   write(*,'(5x,15i4)') (nscl(j,i),j=1,nccl(i))
!end do

!Assign SYVA variables to returned variables
nclass=nseq
classnatm(1:nclass)=nccl(1:nclass)
do i=1,nclass
    classidx(:nccl(i),i)=nscl(:nccl(i),i)
end do

!Determine point group and framework group
call syva_point_group(ng,ni,nsg,ncr,nsr,np,pglabel,nout)
end subroutine



!!-------- Determine HOMO index for single-determinant wavefunction using the safest way. idxHOMO and idxHOMOb are global variables
!Can be used for the case of only GTF information and that containing basis function information
!idxHOMO: HOMO of RHF, or highest SOMO of ROHF, or alpha-HOMO of UHF
!idxHOMOb: beta-HOMO. For RHF, it is equivalent to idxHOMO. If there is no beta electron, idxHOMOb will return 0
!Note that for UHF, actual HOMO energy should be max(MOene(idxHOMO),MOene(idxHOMOb)), because beta-HOMO may be higher than alpha-HOMO
subroutine getHOMOidx
use defvar
if (wfntype==0) then
    do idxHOMO=nmo,1,-1
	    if (nint(MOocc(idxHOMO))==2) exit
    end do
    idxHOMOb=idxHOMO
else if (wfntype==1) then !U
	if (any(MOene/=0)) then !The following way is in principle the safest way, allowing alpha and beta orbitals occur alternately
		ealow=-1E20
		eblow=-1E20
		idxHOMO=0
		idxHOMOb=0
		do imo=1,nmo
			if (nint(MOocc(imo))==0) cycle
			if (MOtype(imo)==1) then
				if (MOene(imo)>ealow) then
					idxHOMO=imo
					ealow=MOene(imo)
				end if
			else if (MOtype(imo)==2) then
				if (MOene(imo)>eblow) then
					idxHOMOb=imo
					eblow=MOene(imo)
				end if
			end if
		end do
    else !If OT is used in CP2K, no orbital energy is available, we have to use the following way and assume all alpha MOs occur prior to beta MOs
		do itmp=nmo,1,-1 !Find the last alpha MO
			if (MOtype(itmp)==1) exit
		end do
		do idxHOMO=itmp,1,-1
			if (nint(MOocc(idxHOMO))==1) exit
		end do
        if (itmp==nmo) then !No beta orbital
			idxHOMOb=0
        else
			do idxHOMOb=nmo,itmp+1,-1
				if (nint(MOocc(idxHOMOb))==1) exit
			end do
        end if
    end if
else if (wfntype==2) then !RO
    do idxHOMO=nmo,1,-1
	    if (nint(MOocc(idxHOMO))==1) exit
    end do
    do idxHOMOb=nmo,1,-1
	    if (nint(MOocc(idxHOMOb))==2) exit
    end do
else
    !write(*,"(a)") " Note: Unable to determine HOMO index because this is not a single-determinant wavefunction"
end if
end subroutine



!!--------- Define fragment corresponding to global variable "frag1". Used by some orbital composition functions
subroutine definefragment
use defvar
use util
character c2000tmp*2000
if (allocated(frag1)) then
	write(*,*) "Atoms in current fragment:"
	write(*,"(13i6)") frag1
	write(*,"(a)") " Input 0 to keep unchanged, or redefine fragment, &
    e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment"
else
	write(*,"(a)") " Input atomic indices to define fragment. &
    e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment"
end if
read(*,"(a)") c2000tmp
if (c2000tmp(1:1)=='0') then
	continue
else if (c2000tmp(1:1)==" ") then
	deallocate(frag1)
else
	if (allocated(frag1)) deallocate(frag1)
	call str2arr(c2000tmp,nfrag1)
	allocate(frag1(nfrag1))
	call str2arr(c2000tmp,nfrag1,frag1)
end if
end subroutine



!!--------- Show very basic information for a range of orbitals (from ibeg to iend)
subroutine showorbinfo(ibeg,iend)
use defvar
integer ibeg,iend,i
do i=1,nmo
	write(*,"(' Orbital:',i5,' Energy(a.u.):',f14.8,' Occ:',f14.8,' Type: ',a)") i,MOene(i),MOocc(i),orbtypename(MOtype(i))
end do
end subroutine



!!--------- Initialize LIBRETA for present wavefunction if haven't
!info=1: Show some notices
!info=2: Silent
!Usually according to "function ifdoESP" to determine if initialize LIBRETA for a real space function
subroutine doinitlibreta(info)
use defvar
use libreta
integer info
if (ifPBC>0) then
	write(*,*) "Error: Currently evaluation of ESP does not support periodic system"
    write(*,*) "Press ENTER button to continue (then Multiwfn may crash)"
    read(*,*)
    pause
end if
if (iuserfunc>=61.and.iuserfunc<=67) return !In this case, the function involves derivative of ESP, while I found the numerical derivate of ESP produced by libreta has noise
if (if_initlibreta==0) then
	if (info==1) then
		write(*,*)
		if (iESPcode==2) write(*,*) "Initializing LIBRETA library (fast version) for ESP evaluation ..."
		if (iESPcode==3) write(*,*) "Initializing LIBRETA library (slow version) for ESP evaluation ..."
		if (nprims>7000) then
			write(*,"(a)") " Note: Number of GTFs of present system is large, if then Multiwfn crashes due to insufficient memory, &
			please change ""iESPcode"" in settings.ini to 1 to use slower ESP evaluation code instead. Alternatively, use a computer with larger memory!"
		end if
    end if
    if (iESPcode==2) then
	    call initlibreta
    else if (iESPcode==3) then
		call initlibreta_slow
    end if
    if_initlibreta=1 !Global variable
    if (info==1) write(*,*) "LIBRETA library has been successfully initialized!"
end if
if (info==1) then
	write(*,"(/,a)") " NOTE: The ESP evaluation code based on LIBRETA library is being used. &
	Please not only cite Multiwfn original paper in your work, but also cite the paper describing the efficient ESP evaluation algorithm adopted by Multiwfn:"
	write(*,"(a)") " Jun Zhang, Tian Lu, Phys. Chem. Chem. Phys., 23, 20323 (2021) DOI: 10.1039/D1CP02805G"
	if (isys==1.and.nthreads>10) then
		write(*,"(a,/)") " Warning!!! In Windows system, it is found that the performance of ESP evaluation code may &
		severely degrade when more than 10 CPU cores are used, therefore 10 cores are used in the following ESP calculation. &
		If you want to pursue better performance by utilizing more cores, please use Linux version instead!"
	end if
end if
end subroutine



!!--------- Initialize for some special functions
subroutine init_func
use defvar
!Local Hartree-Fock exchange energy involves Coulomb matrix generated by libreta and density matrix, so initialize
if (iuserfunc==999.and.allocated(b)) then
    write(*,"(/,a)") " Local Hartree-Fock exchange energy is chosen as user-defined function, &
    it requires Coulomb integral between GTFs, so Libreta library is initialized now..."
    if (iESPcode/=3) then
		write(*,*) "iESPcode in settings.ini is not 3, now is forced to switched to 3"
		iESPcode=3
    end if
    call doinitlibreta(0)
    write(*,"(/,a)") " Note: If local Hartree-Fock exchange energy will be studied using Multiwfn, &
    please note only cite original paper of Multiwfn but also cite Libreta library: J. Chem. Theory Comput., 14, 572 (2018)"
    write(*,*)
    write(*,*) "Generating density matrix between GTFs..."
    call genPprim
end if
end subroutine



!!------------- Convert matrix loaded from ORCA output file to the basis function order in Multiwfn
!Order of ORCA matrices:
!s    
!pz,px,py   
!dz2,dxz,dyz,dx2y2,dxy  
!f0,f+1,f-1,f+2,f-2,f+3,f-3  
!g0,g+1,g-1,g+2,g-2,g+3,g-3,g+4,g-4
!Clearly, we only need to alter the sequence of P shell to make the order as px,py,pz
!  Notice that F(+3) and F(-3) are normalized to -1, therefore the relevant matrix elements must invert sign, &
!while the elements between these two basis functions should keep unchanged due to cancellation
!Simiarly for G(+3,-3,+4,-4) and H(+3,-3,+4,-4)
subroutine ORCAmatconv(ndim,mat)
use defvar
implicit real*8 (a-h,o-z)
real*8 mat(ndim,ndim),tmparr(ndim)

ibas=0
do ish=1,nshell
    ishtype=shtype(ish)
    nshbas=shtype2nbas(ishtype)
    if (ishtype==1) then
        !Reorder row
        tmparr(:)=mat(ibas+1,:) !Backup pz
        mat(ibas+1,:)=mat(ibas+2,:)
        mat(ibas+2,:)=mat(ibas+3,:)
        mat(ibas+3,:)=tmparr(:)
        !Reorder column
        tmparr(:)=mat(:,ibas+1) !Backup pz
        mat(:,ibas+1)=mat(:,ibas+2)
        mat(:,ibas+2)=mat(:,ibas+3)
        mat(:,ibas+3)=tmparr(:)
    else if (ishtype==-3) then !F, invert F(+3,-3)
        mat(ibas+6:ibas+7,:)=-mat(ibas+6:ibas+7,:)
        mat(:,ibas+6:ibas+7)=-mat(:,ibas+6:ibas+7)
    else if (ishtype==-4) then !G, invert G(+3,-3,+4,-4)
        mat(ibas+6:ibas+9,:)=-mat(ibas+6:ibas+9,:)
        mat(:,ibas+6:ibas+9)=-mat(:,ibas+6:ibas+9)
    else if (ishtype==-5) then !H, invert H(+3,-3,+4,-4)
        mat(ibas+6:ibas+9,:)=-mat(ibas+6:ibas+9,:)
        mat(:,ibas+6:ibas+9)=-mat(:,ibas+6:ibas+9)
    end if
    ibas=ibas+nshbas
end do
end subroutine



!!--------- Generate atomic overlap matrix of basis functions for every atom
subroutine genAOMbas(AOMbas)
use defvar
use function
use util
implicit real*8 (a-h,o-z)
real*8 AOMbas(nbasis,nbasis,ncenter),AOMtmp(nbasis,nbasis),basval(nbasis)
real*8 atmspcweight(radpot*sphpot)
type(content) gridatm(radpot*sphpot),gridatmorg(radpot*sphpot)

if (.not.allocated(CObasa_org)) allocate(CObasa_org(nbasis,nbasis))
CObasa_org=CObasa
CObasa=0
do ibas=1,nbasis
    CObasa(ibas,ibas)=1
end do
call CObas2CO(1)
if (.not.allocated(COtr)) allocate(COtr(nprims,nmo))
COtr=transpose(CO) !Global matrix, which will be used in calcbasval for faster calculation

!Decreasing grid quality was found to hinder localization convergence, so do not use this trick
!if (iautointgrid==1) then
!	radpot=45
!	sphpot=170
!end if

!Using grid distance cutoff may hinder localization convergence for e.g. Li6 cluster (deviation of normalization to 1 of this case is quite large)
!So, force to disable this trick. This is probably the contribution of basis function overlap far from atomic center is nonnegligible for sparse system
radcut_old=radcut
radcut=0
call gen1cintgrid(gridatmorg,iradcut) !Generate integration grid

write(*,"(' Radial points:',i5,'    Angular points:',i5,'   Total:',i10,' per center')") radpot,sphpot,radpot*sphpot
call walltime(nwalltime1)

AOMbas=0
ifinish=0
call showprog(0,ncenter)
!Cycle each atom
do iatm=1,ncenter
	gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule
	gridatm%y=gridatmorg%y+a(iatm)%y
	gridatm%z=gridatmorg%z+a(iatm)%z
    gridatm%value=gridatmorg%value
	call gen1cbeckewei(iatm,iradcut,gridatm,atmspcweight)
    
    !$OMP parallel shared(AOMbas) private(ipt,ibas,jbas,AOMtmp,basval,weitmp,weitmp2) num_threads(nthreads)
    AOMtmp=0D0
    !$OMP do schedule(dynamic)
    do ipt=1+iradcut*sphpot,radpot*sphpot
        weitmp=atmspcweight(ipt)*gridatm(ipt)%value
	    !call orbderv(1,1,nbasis,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z,basval)
        call calcbasval(gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z,basval) !Faster than using above line
        !Original version, slower
     !   do ibas=1,nbasis
     !       weitmp2=basval(ibas)*weitmp
		   ! do jbas=ibas,nbasis
			  !  AOMtmp(ibas,jbas)=AOMtmp(ibas,jbas)+basval(jbas)*weitmp2
		   ! end do
	    !end do
        
		do jbas=1,nbasis
            AOMtmp(jbas:nbasis,jbas)=AOMtmp(jbas:nbasis,jbas)+basval(jbas:nbasis)*basval(jbas)*weitmp
	    end do
    end do
    !$OMP end do
    !$OMP CRITICAL
	    AOMbas(:,:,iatm)=AOMbas(:,:,iatm)+AOMtmp(:,:)
    !$OMP end CRITICAL
    !$OMP end parallel
    
    ifinish=ifinish+1
    call showprog(ifinish,ncenter)
end do !End cycling atoms
    
!Original version, slower
!do ibas=1,nbasis
!    do jbas=ibas+1,nbasis
!        AOMbas(jbas,ibas,:)=AOMbas(ibas,jbas,:)
!    end do
!end do
do jbas=1,nbasis
    do ibas=jbas,nbasis
		AOMbas(jbas,ibas,:)=AOMbas(ibas,jbas,:)
	end do
end do

call walltime(nwalltime2)
write(*,"(' Generation of atomic overlap matrix took up',i8,' seconds wall clock time')") nwalltime2-nwalltime1

!Check quality of AOMbas
devmax=0
do ibas=1,nbasis
    tmp=sum(AOMbas(ibas,ibas,:))
    !write(*,"(' Basis function:',i6,'    Normalization:',f16.10)") ibas,tmp
    if (abs(1-tmp)>devmax) devmax=abs(1-tmp)
end do
write(*,"(a,f12.8)") " Maximal deviation of normalization of basis function to unity:",devmax

deallocate(COtr)
CObasa=CObasa_org
CO=CO_org
radcut=radcut_old
end subroutine




!!---- Return electric dipole moment vector, wavefunction information must be available in either b or CObasa
subroutine get_dipole_moment(vec)
use defvar
implicit real*8 (a-h,o-z)
real*8 vec(3)
if (allocated(CObasa)) then
    write(*,*) "Calculating electric dipole moment integral matrix..."
    call genDbas_curr
    xeledip=sum(Dbas(1,:,:)*Ptot(:,:))
    yeledip=sum(Dbas(2,:,:)*Ptot(:,:))
    zeledip=sum(Dbas(3,:,:)*Ptot(:,:))
else if (allocated(b)) then
    write(*,*) "Calculating density matrix based on GTFs..."
    call genPprim
    write(*,*) "Calculating electric dipole moment integral matrix..."
    call genDprim
    xeledip=sum(Dprim(1,:,:)*Ptot_prim(:,:))
    yeledip=sum(Dprim(2,:,:)*Ptot_prim(:,:))
    zeledip=sum(Dprim(3,:,:)*Ptot_prim(:,:))
end if
xnucdip=0
ynucdip=0
znucdip=0
do iatm=1,ncenter
    xnucdip=xnucdip+a(iatm)%x*a(iatm)%charge
    ynucdip=ynucdip+a(iatm)%y*a(iatm)%charge
    znucdip=znucdip+a(iatm)%z*a(iatm)%charge
end do
vec(1)=xeledip+xnucdip
vec(2)=yeledip+ynucdip
vec(3)=zeledip+znucdip
end subroutine


!!------ Generate "fragatm" array containing all atoms in the present system
subroutine genfragatm
use defvar
if (allocated(fragatm)) deallocate(fragatm)
allocate(fragatm(ncenter))
nfragatm=ncenter
forall (i=1:nfragatm) fragatm(i)=i
ifragcontri=0 !Fragment has not been manually defined by user
end subroutine


!!-------- Export all internal geometry parameters to int_coord.txt in current folder
subroutine showgeomparam(outname)
use defvar
use util
character(len=*) outname
integer nneigh(ncenter) !Number of neighbours (the number of atoms connect to this atom)
integer neigh(maxneigh,ncenter) !neigh(1:nneigh(i),i) is list of neighbouring atom indices of atom i
integer nbond !Number of bonds
integer,allocatable :: bond(:,:) !bond(1/2,i) are indices of the two atoms involved in bond i
integer nangle !Number of angles
integer,allocatable :: angle(:,:) !angle(1/2/3,i) are indices of the three atoms involved in angle i
integer ndih !Number of dihedrals
integer,allocatable :: dih(:,:) !dih(1/2/3/4,i) are indices of the four atoms involved in dihedral i

open(10,file=outname,status="replace")

!Determining neighbours
nneigh=0
do iatm=1,ncenter
    do jatm=1,ncenter
        if (jatm==iatm) cycle
        if (distmat(iatm,jatm)<( covr(a(iatm)%index)+covr(a(jatm)%index) )*bondcrit) then
            nneigh(iatm)=nneigh(iatm)+1
            neigh(nneigh(iatm),iatm)=jatm
        end if
    end do
end do

!Determining bond list
do itime=1,2
    nbond=0
    do iatm=1,ncenter
        do jatm=iatm+1,ncenter
            if (any(neigh(1:nneigh(iatm),iatm)==jatm)) then
                nbond=nbond+1
                if (itime==2) then
                    bond(1,nbond)=iatm
                    bond(2,nbond)=jatm
                end if
            end if
        end do
    end do
    if (itime==1) allocate(bond(2,nbond))
end do
write(10,"(' Number of bonds:',i8)") nbond
do ibond=1,nbond
    iatm=bond(1,ibond)
    jatm=bond(2,ibond)
    write(10,"(' #',i5,'      Atoms:',2i5,'    Distance:',f12.6,' Angstrom')") ibond,iatm,jatm,distmat(iatm,jatm)*b2a
end do

!Determining angle list
do itime=1,2
    nangle=0
    do icen=1,ncenter
        if (nneigh(icen)<2) cycle
        do idx=1,nneigh(icen)
            iatm=neigh(idx,icen)
            do jdx=idx+1,nneigh(icen)
                jatm=neigh(jdx,icen)
                nangle=nangle+1
                if (itime==2) then
                    angle(1,nangle)=iatm
                    angle(2,nangle)=icen
                    angle(3,nangle)=jatm
                end if
            end do
        end do
    end do
    if (itime==1) allocate(angle(3,nangle))
end do
!Sort angle array, making index from small to large, from left to right
call sortidxlist(angle,3,nangle)
write(10,*)
write(10,"(' Number of angles:',i8)") nangle
do iangle=1,nangle
    iatm=angle(1,iangle)
    jatm=angle(2,iangle)
    katm=angle(3,iangle)
    angval=xyz2angle(a(iatm)%x,a(iatm)%y,a(iatm)%z,a(jatm)%x,a(jatm)%y,a(jatm)%z,a(katm)%x,a(katm)%y,a(katm)%z)
    write(10,"(' #',i5,'      Atoms:',3i5,'    Angle:',f12.6,' degree')") iangle,iatm,jatm,katm,angval
end do

!Determining dihedral list. iatm-icen-jcen-jatm
nlinear=0
do itime=1,2
    ndih=0
    do ibond=1,nbond
        icen=bond(1,ibond)
        jcen=bond(2,ibond)
        neiA=nneigh(icen)
        neiB=nneigh(jcen)
        if (neiA==0.or.neiB==0) cycle
        do idx=1,neiA
            iatm=neigh(idx,icen)
            if (iatm==jcen) cycle
            do jdx=1,neiB
                jatm=neigh(jdx,jcen)
                if (jatm==icen) cycle
                if (iatm==jatm) cycle
                angdev1=180-atomang(iatm,icen,jcen)
                angdev2=180-atomang(icen,jcen,jatm)
                if (angdev1<1.or.angdev2<1) then
                    nlinear=nlinear+1
                    cycle
                end if
                ndih=ndih+1
                if (itime==2) then
                    if (jatm>iatm) then
                        dih(1,ndih)=iatm
                        dih(2,ndih)=icen
                        dih(3,ndih)=jcen
                        dih(4,ndih)=jatm
                    else !Require the index of the fourth atom is larger than that of the first atom
                        dih(1,ndih)=jatm
                        dih(2,ndih)=jcen
                        dih(3,ndih)=icen
                        dih(4,ndih)=iatm
                    end if
                end if
            end do
        end do
    end do
    if (itime==1) allocate(dih(4,ndih))
end do
if (nlinear>0) then
    write(*,"(a,i5,a)") " Note:",nlinear," dihedrals deviate from linear less than 1 degree, they are ignored"
end if

!Sort dih array, making index from small to large, from left to right
call sortidxlist(dih,4,ndih)
write(10,*)
write(10,"(' Number of dihedrals:',i8)") ndih
do idih=1,ndih
    iatm=dih(1,idih)
    jatm=dih(2,idih)
    katm=dih(3,idih)
    latm=dih(4,idih)
    dihval=xyz2dih(a(iatm)%x,a(iatm)%y,a(iatm)%z,a(jatm)%x,a(jatm)%y,a(jatm)%z,a(katm)%x,a(katm)%y,a(katm)%z,a(latm)%x,a(latm)%y,a(latm)%z)
    write(10,"(' #',i5,'      Atoms:',4i5,'    Dihedral:',f12.6,' degree')") idih,iatm,jatm,katm,latm,dihval
end do

close(10)

end subroutine




!!------ Generate Z-matrix. The each line of Zmat corresponds to an entry of the matrix, atomic indices are recorded
!If ierror=1, then the generation is failed (e.g. C2H2). =0 means normal
subroutine genZmat(Zmat,ierror)
use defvar
use util
integer Zmat(ncenter,3),ierror
integer nneigh(ncenter) !Number of neighbours (the number of atoms connected to this atom)
integer neigh(maxneigh,ncenter) !neigh(1:nneigh(i),i) is list of neighbouring atom indices of atom i

call genneighlist(1,1,neigh,nneigh) !Determining neighbours according to connmat

ierror=0
Zmat=0 !Element of 0 means unavailable

!Content of a line: iatm i1 i2 i3
!Will determine suitable i1, then i2, then i3
!When determining the suitable one w.r.t. current atom p, first search an atom with index less than iatm and connects to p, &
!if nothing is found, use an atom with index less than iatm and nearest to p
!Three atoms should not lie nearly a straight line, the angle must be >0.5 and <179.5 degree

do iatm=2,ncenter
    do jatm=1,iatm-1 !If an atom has index less than iatm, use it
        if (any(neigh(1:nneigh(iatm),iatm)==jatm)) then
            i1=jatm
            exit
        end if
    end do
    if (jatm==iatm) then !No atom connected to iatm, find the atom nearest to the current one
        distmin=1E10
        do jatm=1,iatm-1
            dist=atomdist(iatm,jatm)
            if (dist<distmin) then
                i1=jatm
                distmin=dist
            end if
        end do
    end if
    Zmat(iatm,1)=i1
    
    if (iatm>=3) then
        do jatm=1,iatm-1 !If an atom has index less than i1, use it
            if (jatm==i1) cycle
            if (any(neigh(1:nneigh(i1),i1)==jatm)) then
				if (ncenter>3) then
					tmpang=xyz2angle(a(iatm)%x,a(iatm)%y,a(iatm)%z,a(i1)%x,a(i1)%y,a(i1)%z,a(jatm)%x,a(jatm)%y,a(jatm)%z)
					if (tmpang<0.5D0.or.tmpang>179.5D0) cycle
                end if
                i2=jatm
                exit
            end if
        end do
        if (jatm==iatm) then !No atom connected to iatm, find the atom nearest to the current one
            distmin=1E10
            do jatm=1,iatm-1
                if (jatm==i1) cycle
                dist=atomdist(i1,jatm)
                if (dist<distmin) then
					if (ncenter>3) then
						tmpang=xyz2angle(a(iatm)%x,a(iatm)%y,a(iatm)%z,a(i1)%x,a(i1)%y,a(i1)%z,a(jatm)%x,a(jatm)%y,a(jatm)%z)
						if (tmpang<0.5D0.or.tmpang>179.5D0) cycle
                    end if
                    i2=jatm
                    distmin=dist
                end if
            end do
            if (distmin==1E10) then !distmin was unchanged, no atom satisfies condition
				ierror=1
                return
            end if
        end if
        Zmat(iatm,2)=i2
    
        if (iatm>=4) then
            do jatm=1,iatm-1 !If an atom has index less than iatm, use it
                if (jatm==i1.or.jatm==i2) cycle
                if (any(neigh(1:nneigh(i2),i2)==jatm)) then
					tmpang=xyz2angle(a(i1)%x,a(i1)%y,a(i1)%z,a(i2)%x,a(i2)%y,a(i2)%z,a(jatm)%x,a(jatm)%y,a(jatm)%z)
					if (tmpang<0.5D0.or.tmpang>179.5D0) cycle
					i3=jatm
					exit
                end if
            end do
            if (jatm==iatm) then !No atom connected to iatm, find the atom nearest to the current one
                distmin=1E10
                do jatm=1,iatm-1
                    if (jatm==i1.or.jatm==i2) cycle
                    dist=atomdist(i2,jatm)
                    if (dist<distmin) then
						tmpang=xyz2angle(a(i1)%x,a(i1)%y,a(i1)%z,a(i2)%x,a(i2)%y,a(i2)%z,a(jatm)%x,a(jatm)%y,a(jatm)%z)
						if (tmpang<0.5D0.or.tmpang>179.5D0) cycle
                        i3=jatm
                        distmin=dist
                    end if
                end do
				if (distmin==1E10) then
					ierror=1
					return
				end if
            end if
            Zmat(iatm,3)=i3
        end if
    end if
	!write(*,*) iatm,Zmat(iatm,1:3)
end do
end subroutine




!!---- Print information to let user to input orbital label
!Only supports single-determinant wavefunction, except for RO, and basis function should be available
subroutine orblabsel_prompt
use defvar
if ((wfntype==0.or.wfntype==1).and.allocated(CObasa)) then
	if (wfntype==0) then
		write(*,*) "You can also input for example: (quotation marks should not be inputted)"
		write(*,*) """h"" for HOMO, ""h-5"" for HOMO-5, ""l"" for LUMO, ""l+3"" for LUMO+3"
	else
		write(*,*) "You can also input for example: (quotation marks should not be inputted)"
		write(*,*) """ha"" for HOMO(alpha), ""hb"" for HOMO(beta), ""hb-5"" for HOMO(beta)-5"
		write(*,*) """la"" for LUMO(alpha), ""lb"" for LUMO(beta), ""la+3"" for LUMO(alpha)+3"
	end if
end if
end subroutine

!!------- Return orbital index according to inputted orbital label
!Only supports single-determinant wavefunction, except for RO, and basis function should be available
!iorb is orbital index. If failed to determine, return 0
!Suggested usage:
!if (index(c200tmp,'h')==0.and.index(c200tmp,'l')==0) then 
!   ... directly load orbital index
!else
!   call orblabsel(c200tmp,iorb)
!	if (ishowmo==0) then
!		write(*,*) "Error: The orbital label you inputted is wrong! Please double check"
!		...cycle
!	end if
!end if
subroutine orblabsel(orblab,iorb)
use defvar
character(len=*) orblab
integer iorb,tmpval
if (index(orblab,',')/=0) then
    write(*,*) "Error: In this selection mode, you can select only one orbital each time"
    iorb=0
    return
end if
call getHOMOidx
tmpval=0
orblab=adjustl(orblab)
if (orblab(1:2)=="ha") then
    if (orblab(3:3)=='-') read(orblab(3:),*) tmpval
    iorb=idxHOMO+tmpval
else if (orblab(1:2)=="hb") then
    if (orblab(3:3)=='-') read(orblab(3:),*) tmpval
    iorb=idxHOMOb+tmpval
else if (orblab(1:2)=="h ".or.orblab(1:2)=="h-") then
    if (orblab(2:2)=='-') read(orblab(2:),*) tmpval
    iorb=idxHOMO+tmpval
else if (orblab(1:2)=="la") then
    if (orblab(3:3)=='+') read(orblab(3:),*) tmpval
    iorb=idxHOMO+1+tmpval
else if (orblab(1:2)=="lb") then
    if (orblab(3:3)=='+') read(orblab(3:),*) tmpval
    iorb=idxHOMOb+1+tmpval
else if (orblab(1:2)=="l ".or.orblab(1:2)=="l+") then
    if (orblab(2:2)=='+') read(orblab(2:),*) tmpval
    iorb=idxHOMO+1+tmpval
else
    iorb=0
end if
if (iorb<1.or.iorb>nmo) iorb=0
end subroutine




!!----------- Transform MO energies to Fock matrix
!The numerical accuracy may be marginally lower than that directly loaded from e.g. 47
!istatus is returned variable. It is 0 by default, when one or more orbital energies are zero and thus cannot generate Fock, istatus=1
subroutine MOene2Fmat(istatus)
use defvar
use util
implicit real*8 (a-h,o-z)
real*8 tmpmat(nbasis,nbasis),Emat(nbasis,nbasis),Cinv(nbasis,nbasis)
integer istatus

istatus=0
if (.not.allocated(CObasa)) then
	write(*,*) "Error: Basis function information is not available!"
    write(*,*) "Press ENTER button to return"
    read(*,*)
    istatus=1
    return
end if
nzero=count(MOene==0)
if (nzero>0) then
    if (nzero==nbasis) then
        write(*,"(a)") " Error: All orbital energies are zero, hence Fock/KS matrix &
        cannot be generated based on orbital energies and expansion coefficients. Commonly this is because the orbital energies were not &
        solved at all during your quantum chemistry calculation. For example, the current orbitals are not molecular orbitals, or you enabled OT in CP2K"
    else
        write(*,"(a,i6,a)") " Error:",nzero," orbital energies are zero, hence Fock/KS matrix &
        cannot be generated based on orbital energies and expansion coefficients. Commonly this is because there is linear dependency problem &
        in your quantum chemistry calculation, and thus some basis functions are automatically removed by your quantum chemistry code. &
        To tackle this problem, removing diffuse functions if they have been used; if you are a Gaussian user, also add IOp(3/32=2) to fully avoid &
        eliminating linear-dependent basis functions"
    end if
    write(*,*) "Press ENTER button to return"
    read(*,*)
    istatus=1
    return
end if

write(*,*) "Generating Fock/KS matrix based on orbital energies and LCAO coefficients ..."

!Because FC=SCE, thus F=SCE(C)^-1
if (allocated(FmatA)) deallocate(FmatA)
allocate(FmatA(nbasis,nbasis))
Cinv=invmat(CObasa,nbasis)
Emat=0
forall(i=1:nbasis) Emat(i,i)=MOene(i)
!FmatA=matmul(matmul(matmul(Sbas,CObasa),Emat),Cinv) !Slower code
tmpmat=matmul_blas(Sbas,CObasa,nbasis,nbasis,0,0) !Parallel MKL
tmpmat=matmul_blas(tmpmat,Emat,nbasis,nbasis,0,0)
FmatA=matmul_blas(tmpmat,Cinv,nbasis,nbasis,0,0)

if (wfntype==1.or.wfntype==4) then !Beta part
    if (allocated(FmatB)) deallocate(FmatB)
    allocate(FmatB(nbasis,nbasis))
    Cinv=invmat(CObasb,nbasis)
    Emat=0
    forall(i=1:nbasis) Emat(i,i)=MOene(i+nbasis)
    tmpmat=matmul_blas(Sbas,CObasb,nbasis,nbasis,0,0) !Parallel MKL
    tmpmat=matmul_blas(tmpmat,Emat,nbasis,nbasis,0,0)
    FmatB=matmul_blas(tmpmat,Cinv,nbasis,nbasis,0,0)
end if

write(*,*) "Fock/KS matrix has been successfully generated!"

!Compare with that loaded from .47
!tmpmat=FmatA
!call loadFockfile(istatus)
!!call showmatgau(FmatA,form="f14.8")
!write(*,"(' Maximum deviation to loaded Fock matrix',f20.10)") maxval(abs(tmpmat-FmatA))
end subroutine



!!----------- Find element index from element name. If idx returns 0, means the index was not found
subroutine elename2idx(name,idx)
use defvar
use util
character(len=*) name
integer idx
idx=0
call lc2uc(name(1:1)) !Convert to upper case
call uc2lc(name(2:2)) !Convert to lower case
do iele=1,nelesupp
	if ( name(1:2)==ind2name(iele) ) then
		idx=iele
		return
	end if
end do
end subroutine