diff --git a/diagonalization/chani.F90 b/diagonalization/chani.F90 deleted file mode 100644 index dc814df2afe564e7d5a106e445f948f66585f8bc..0000000000000000000000000000000000000000 --- a/diagonalization/chani.F90 +++ /dev/null @@ -1,441 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -MODULE m_chani -CONTAINS - SUBROUTINE chani(M,Neigd, Myid,Np,Sub_comm,mpi_comm, & - Eig,Num,hamOvlp,zMat) - ! - !---------------------------------------------------- - !- Parallel eigensystem solver - driver routine; gb99 - ! Uses the SCALAPACK for the actual diagonalization - ! - ! m ........ actual (=leading) dimension of full a & b matrices - ! must be problem size, as input a, b are one-dimensional - ! and shall be redistributed to two-dimensional matrices - ! actual (=leading) dimension of eigenvector z(,) - ! n ........ number of columns of full (sub)matrix ( about n/np) - ! neigd..... second (=column) dimension of eigenvector matrix - ! myid ..... id of node (normally irank) - ! np ....... number of processors (normally isize) - ! SUB_COMM.. communicator for MPI - ! a,b .... packed (sub)matrices, here expanded to non-packed - ! z,eig .... eigenvectors and values, output - ! num ...... number of ev's searched (and found) on this node - ! On input, overall number of ev's searched, - ! On output, local number of ev's found - ! - !---------------------------------------------------- - ! - !#include"cpp_arch.h" -#include"cpp_double.h" - USE m_juDFT - USE m_types - USE m_subredist1 - USE m_subredist2 - IMPLICIT NONE - INCLUDE 'mpif.h' - - INTEGER, INTENT (IN) :: neigd,m - INTEGER, INTENT (IN) :: SUB_COMM,np,myid,mpi_comm - INTEGER, INTENT (INOUT) :: num - REAL, INTENT (OUT) :: eig(neigd) - TYPE(t_hamOvlp),INTENT(INOUT) :: hamOvlp - TYPE(t_zMat),INTENT(INOUT) :: zMat - - !... Local variables - ! - INTEGER nc,ic,ir,n_sym,jsym,num_j,icm,n_bound - INTEGER i ,j ,k ,l, info, i_plus, i_dim - INTEGER nprow,npcol,myrowssca, nsq, nprow2, npcol2 - INTEGER myrow,mycol,mycolssca, ierr, me,ierr2 - INTEGER iamblacs,myrowblacs,mycolblacs,npblacs - INTEGER ictxtblacs,err - INTEGER, ALLOCATABLE :: iwork(:),iusermap(:,:) - INTEGER, ALLOCATABLE :: iblacsnums(:),ihelp(:) - REAL, ALLOCATABLE :: dwork(:) - REAL, ALLOCATABLE :: eigvec_r(:,:) - COMPLEX, ALLOCATABLE :: eigvec_c(:,:) - REAL, ALLOCATABLE :: rwork(:) - INTEGER :: lrwork - - ! - ! ScaLAPACK things - CHARACTER (len=1) :: uplo - INTEGER, PARAMETER :: dlen_=9 - INTEGER :: desca(dlen_),desceigv(dlen_) - INTEGER :: num1,num2,liwork,lwork2,np0,mq0 - INTEGER :: iceil, numroc, nn, nb - INTEGER, ALLOCATABLE :: ifail(:), iclustr(:) - REAL :: abstol,orfac=1.E-4,CPP_LAPACK_slamch - REAL,ALLOCATABLE :: eig2(:), gap(:) - REAL, ALLOCATABLE :: asca_r(:,:), bsca_r(:,:),work2_r(:) - COMPLEX, ALLOCATABLE :: asca_c(:,:), bsca_c(:,:),work2_c(:) - - EXTERNAL iceil, numroc - EXTERNAL CPP_LAPACK_slamch, descinit - EXTERNAL blacs_pinfo, blacs_gridinit - EXTERNAL MPI_COMM_DUP - - ! - ! determine actual number of columns of input matrices A and B - ! nc is number of columns the local processor will get, must be <=n - ! - nc = 0 - DO i = myid+1, m, np - nc = nc + 1 - ENDDO - !IF (nc.GT.n) THEN - ! WRITE (6,*) myid,'gets more columns than allowed:' - ! WRITE (6,*) myid,'will get',nc,' columns, only',n,' allowed' - ! CALL juDFT_error("chani: nc > n",calledby ="chani") - !ENDIF - ! - ! determine block size - ! - nb = 20 - IF (m.GT.2048) nb = 30 !2 - IF (m.GT.8*2048) nb = 60 !4 - ! - ! transfer to ScaLAPACK block cyclic 2d distribution with - ! block-size nb - ! - ! compute processor grid, as square as possible - ! If not square with more rows than columns - nsq = INT(SQRT(REAL(np))) - gridloop:DO j=nsq,1,-1 - IF ( (np/j) * j == np) THEN - npcol = np/j - nprow = j - EXIT gridloop - ENDIF - ENDDO gridloop - ! - ! An nprow*npcol processor grid will be created - ! Row and column index myrow, mycol of this processor in the grid - ! and distribution of A and B in ScaLAPACK - ! The local processor will get myrowssca rows and mycolssca columns - ! of A and B - ! - myrow = myid/npcol ! my row number in the BLACS nprow*npcol grid - mycol = myid -(myid/npcol)*npcol ! my column number in the BLACS nprow*npcol grid - ! - ! Now allocate Asca to put the elements of Achi or receivebuffer to - ! - myrowssca=(m-1)/(nb*nprow)*nb+& - MIN(MAX(m-(m-1)/(nb*nprow)*nb*nprow-nb*myrow,0),nb) - ! Number of rows the local process gets in ScaLAPACK distribution - mycolssca=(m-1)/(nb*npcol)*nb+ & - MIN(MAX(m-(m-1)/(nb*npcol)*nb*npcol-nb*mycol,0),nb) - ! Number of columns the local process gets in ScaLAPACK distribution - if (hamOvlp%l_real) THEN - ALLOCATE ( asca_r(myrowssca,mycolssca), stat=err ) - else - ALLOCATE ( asca_c(myrowssca,mycolssca), stat=err ) - endif - IF (err.NE.0) THEN - WRITE (*,*) 'In chani an error occured during the allocation of' - WRITE (*,*) 'asca: ',err,' size: ',myrowssca*mycolssca - CALL juDFT_error("allocte in chani",calledby ="chani") - ENDIF - if (hamOvlp%l_real) THEN - asca_r=0.0 - ALLOCATE ( bsca_r(myrowssca,mycolssca), stat=err ) - else - asca_c=0.0 - ALLOCATE ( bsca_c(myrowssca,mycolssca), stat=err ) - endif - IF (err.NE.0) THEN - WRITE (*,*) 'In chani an error occured during the allocation of' - WRITE (*,*) 'bsca :',err - CALL juDFT_error("allocate in chani",calledby ="chani") - ENDIF - if (hamOvlp%l_real) THEN - bsca_r=0.0 - CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_r=hamovlp%a_r,asca_r=asca_r) - CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_r=hamovlp%b_r,asca_r=bsca_r) - else - bsca_c=0.0 - CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_c=hamovlp%a_c,asca_c=asca_c) - CALL subredist1(m,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ,achi_c=hamovlp%b_c,asca_c=bsca_c) - end if - CALL BLACS_PINFO(iamblacs,npblacs) ! iamblacs = local process rank (e.g. myid) - ! npblacs = number of available processes - IF (npblacs /= np) THEN - !gb WRITE (6,*) 'Number of procs differ between', - !gb & ' BLACS and MPI' - IF (np > npblacs) THEN - WRITE(6,*) 'Have to take SUB_COMM for BLACS grid!' ! this should not happen - ENDIF - ENDIF - ALLOCATE( iblacsnums(np) ) - iblacsnums=-2 - ALLOCATE (ihelp(np)) - ihelp=-2 - ihelp(myid+1)=iamblacs ! Get the Blacs id corresponding to the MPI id - CALL MPI_ALLREDUCE(ihelp, iblacsnums, np,MPI_INTEGER,MPI_MAX,SUB_COMM,ierr) - IF (ierr.NE.0) THEN - WRITE (6,*) 'Error in allreduce for BLACS nums' - CALL juDFT_error('Error in allreduce for BLACS nums',calledby ='chani') - ENDIF - ! iblacsnums(i) is the BLACS-process number of MPI-process i-1 - ALLOCATE ( iusermap(nprow,npcol) ) ! Usermap for BLACS_GRIDMAP - k = 1 - DO i = 1, nprow - DO j = 1, npcol - iusermap(i,j) = iblacsnums(k) - k = k + 1 - ENDDO - ENDDO - !gb WRITE(*,*) np,iblacsnums - ! - ! Take a copy of SUB_COMM as system context for BLACS_GRIDMAP - ! We cannot take SUB_COMM itself as it will be overwritten with the - ! BLACS context for the grid comstructed - ! - ! CALL MPI_COMM_DUP(SUB_COMM,ictxtblacs,ierr) - !I do not understand this, perhaps we should use MPI_COMM_WRLD??? - ! CALL MPI_COMM_DUP(MPI_COMM_WRLD,ictxtblacs,ierr) - CALL MPI_COMM_DUP(MPI_COMM,ictxtblacs,ierr) - IF (ierr /=0 ) THEN - WRITE(6,*) 'MPI_COMM_DUP failed' - CALL juDFT_error('MPI_COMM_DUP failed',calledby='chani') - ENDIF - CALL BLACS_GRIDMAP(ictxtblacs,iusermap,nprow,nprow,npcol) - ! Now control, whether the BLACS grid is the one we wanted - CALL BLACS_GRIDINFO(ictxtblacs, nprow2,npcol2,myrowblacs,mycolblacs) - IF (nprow2 /= nprow) THEN - WRITE(6,*) 'Wrong number of rows in BLACS grid' - WRITE(6,*) 'nprow=',nprow,' nprow2=',nprow2 - ierr = -nprow2 - CALL juDFT_error('Wrong number of rows in BLACS grid',calledby= 'chani') - ENDIF - IF (npcol2 /= npcol) THEN - WRITE(6,*) 'Wrong number of columns in BLACS grid' - WRITE(6,*) 'npcol=',npcol,' npcol2=',npcol2 - ierr = -npcol2 - CALL juDFT_error('Wrong number of columns in BLACS grid', calledby ='chani') - - ENDIF - IF (myrowblacs.NE.myrow) THEN - WRITE(6,*) 'My row in BLACS is',myrowblacs - WRITE(6,*) 'It is not',myrow - !CALL CPP_flush(6) - ENDIF - IF (mycolblacs.NE.mycol) THEN - WRITE(6,*) 'My column in BLACS is',mycolblacs - WRITE(6,*) 'It is not',mycol - !CALL CPP_flush(6) - ENDIF - ! - CALL descinit(DESCA,m,m,nb,nb,0,0,ictxtblacs,myrowssca,ierr) - CALL descinit(DESCeigv,m,m,nb,nb,0,0,ictxtblacs, myrowssca,ierr) - - abstol=2.0*CPP_LAPACK_slamch('S') ! PDLAMCH gave an error on ZAMpano - ALLOCATE ( eig2(m), stat=err ) ! The eigenvalue array for ScaLAPACK - IF (err.NE.0) THEN - WRITE (*,*) 'In chani an error occured during the allocation of' - WRITE (*,*) 'eig2 :',err - CALL juDFT_error('Failed to allocated "eig2"', calledby ='chani') - ENDIF - ! write(*,*) 'c :',myrowssca,mycolssca,desceigv - if (hamovlp%l_real) THEN - ALLOCATE ( eigvec_r(myrowssca,mycolssca), stat=err ) ! Eigenvectors for ScaLAPACK - else - ALLOCATE ( eigvec_c(myrowssca,mycolssca), stat=err ) ! Eigenvectors for ScaLAPACK - end if - IF (err.NE.0) THEN - WRITE (*,*) 'In chani an error occured during the allocation of' - WRITE (*,*) 'eigvec :',err - CALL juDFT_error('Failed to allocated "eigvec"', calledby ='chani') - ENDIF - ! - nn=MAX(MAX(m,nb),2) - np0=numroc(nn,nb,0,0,nprow) - mq0=numroc(MAX(MAX(neigd,nb),2),nb,0,0,npcol) - if (hamovlp%l_real) THEN - lwork2=5*m+MAX(5*nn,np0*mq0+2*nb*nb)+ iceil(neigd,nprow*npcol)*nn - ALLOCATE ( work2_r(lwork2+10*m), stat=err ) ! Allocate more in case of clusters - else - lwork2=m+MAX(nb*(np0+1),3) - ALLOCATE ( work2_c(lwork2), stat=err ) - endif - IF (err.NE.0) THEN - WRITE (*,*) 'work2 :',err,lwork2 - CALL juDFT_error('Failed to allocated "work2"', calledby ='chani') - ENDIF - - liwork=6*MAX(MAX(m,nprow*npcol+1),4) - ALLOCATE ( iwork(liwork), stat=err ) - IF (err.NE.0) THEN - WRITE (*,*) 'iwork :',err,liwork - CALL juDFT_error('Failed to allocated "iwork"', calledby ='chani') - ENDIF - ALLOCATE ( ifail(m), stat=err ) - IF (err.NE.0) THEN - WRITE (*,*) 'ifail :',err,m - CALL juDFT_error('Failed to allocated "ifail"', calledby ='chani') - ENDIF - ALLOCATE ( iclustr(2*nprow*npcol), stat=err ) - IF (err.NE.0) THEN - WRITE (*,*) 'iclustr:',err,2*nprow*npcol - CALL juDFT_error('Failed to allocated "iclustr"', calledby ='chani') - ENDIF - ALLOCATE ( gap(nprow*npcol), stat=err ) - IF (err.NE.0) THEN - WRITE (*,*) 'gap :',err,nprow*npcol - CALL juDFT_error('Failed to allocated "gap"', calledby ='chani') - ENDIF - ! - ! Compute size of workspace - ! - if (hamovlp%l_real) THEN - uplo='U' - CALL CPP_LAPACK_pdsygvx(1,'V','I','U',m,asca_r,1,1,desca,bsca_r,1,1, desca,& - 0.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_r,1,1,& - desceigv,work2_r,-1,iwork,-1,ifail,iclustr, gap,ierr) - IF ( work2_r(1).GT.lwork2) THEN - lwork2 = work2_r(1) - DEALLOCATE (work2_r) - ALLOCATE ( work2_r(lwork2+20*m), stat=err ) ! Allocate even more in case of clusters - IF (err.NE.0) THEN - WRITE (*,*) 'work2 :',err,lwork2 - CALL juDFT_error('Failed to allocated "work2"', calledby ='chani') - ENDIF - ENDIF - else - lrwork=4*m+MAX(5*nn,np0*mq0)+ iceil(neigd,nprow*npcol)*nn - ! Allocate more in case of clusters - ALLOCATE(rwork(lrwork+10*m), stat=ierr) - IF (err /= 0) THEN - WRITE (*,*) 'ERROR: chani.F: Allocating rwork failed' - CALL juDFT_error('Failed to allocated "rwork"', calledby ='chani') - ENDIF - - CALL CPP_LAPACK_pzhegvx(1,'V','I','U',m,asca_c,1,1,desca,bsca_c,1,1, desca,& - 0.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_c,1,1,& - desceigv,work2_c,-1,rwork,-1,iwork,-1,ifail,iclustr,& - gap,ierr) - IF (ABS(work2_c(1)).GT.lwork2) THEN - lwork2=work2_c(1) - DEALLOCATE (work2_c) - ALLOCATE (work2_c(lwork2), stat=err) - IF (err /= 0) THEN - WRITE (*,*) 'ERROR: chani.F: Allocating rwork failed:',lwork2 - CALL juDFT_error('Failed to allocated "work2"', calledby ='chani') - ENDIF - ENDIF - IF (rwork(1).GT.lrwork) THEN - lrwork=rwork(1) - DEALLOCATE(rwork) - ! Allocate even more in case of clusters - ALLOCATE (rwork(lrwork+20*m), stat=err) - IF (err /= 0) THEN - WRITE (*,*) 'ERROR: chani.F: Allocating rwork failed: ', lrwork+20*m - CALL juDFT_error('Failed to allocated "rwork"', calledby ='chani') - ENDIF - ENDIF - endif - IF (iwork(1) .GT. liwork) THEN - liwork = iwork(1) - DEALLOCATE (iwork) - ALLOCATE (iwork(liwork), stat=err) - IF (err /= 0) THEN - WRITE (*,*) 'ERROR: chani.F: Allocating iwork failed: ',liwork - CALL juDFT_error('Failed to allocated "iwork"', calledby ='chani') - ENDIF - ENDIF - ! - ! Now solve generalized eigenvalue problem - ! - CALL timestart("SCALAPACK call") - if (hamovlp%l_real) THEN - CALL CPP_LAPACK_pdsygvx(1,'V','I','U',m,asca_r,1,1,desca,bsca_r,1,1, desca,& - 1.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_r,1,1,& - desceigv,work2_r,lwork2,iwork,liwork,ifail,iclustr,& - gap,ierr) - else - CALL CPP_LAPACK_pzhegvx(1,'V','I','U',m,asca_c,1,1,desca,bsca_c,1,1, desca,& - 1.0,1.0,1,num,abstol,num1,num2,eig2,orfac,eigvec_c,1,1,& - desceigv,work2_c,lwork2,rwork,lrwork,iwork,liwork,& - ifail,iclustr,gap,ierr) - DEALLOCATE(rwork) - endif - CALL timestop("SCALAPACK call") - IF (ierr .NE. 0) THEN - IF (ierr /= 2) WRITE (6,*) myid,' error in pzhegvx/pdsygvx, ierr=',ierr - IF (ierr <= 0) WRITE (6,*) myid,' illegal input argument' - IF (MOD(ierr,2) /= 0) THEN - WRITE (6,*) myid,'some eigenvectors failed to converge' - eigs: DO i = 1, neigd - IF (ifail(i) /= 0) THEN - WRITE (6,*) myid,' eigenvector',ifail(i), 'failed to converge' - ELSE - EXIT eigs - ENDIF - ENDDO eigs - !CALL CPP_flush(6) - ENDIF - IF (MOD(ierr/4,2).NE.0) THEN - WRITE(6,*) myid,' only',num2,' eigenvectors converged' - !CALL CPP_flush(6) - ENDIF - IF (MOD(ierr/8,2).NE.0) THEN - WRITE(6,*) myid,' PDSTEBZ failed to compute eigenvalues' - CALL judft_error("SCALAPACK failed to solve eigenvalue problem",calledby="chani.F90") - ENDIF - IF (MOD(ierr/16,2).NE.0) THEN - WRITE(6,*) myid,' B was not positive definite, Cholesky failed at',ifail(1) - CALL judft_error("SCALAPACK failed: B was not positive definite",calledby="chani.F90") - ENDIF - ENDIF - IF (num2 < num1) THEN - IF (myid ==0) THEN - WRITE(6,*) 'Not all eigenvalues wanted are found' - WRITE(6,*) 'number of eigenvalues/vectors wanted',num1 - WRITE(6,*) 'number of eigenvalues/vectors found',num2 - !CALL CPP_flush(6) - ENDIF - ENDIF - ierr=0 - CALL BLACS_GRIDEXIT(ictxtblacs,ierr) - IF (ierr /= 0 ) THEN - WRITE (6,*) 'blacs_gridexit failed, ierr=',ierr - ENDIF - DEALLOCATE(iwork) - DEALLOCATE(ifail) - DEALLOCATE(iclustr) - DEALLOCATE(gap) - ! - ! Put those eigenvalues expected by chani to eig, i.e. for - ! process i these are eigenvalues i+1, np+i+1, 2*np+i+1... - ! Only num=num2/np eigenvalues per process - ! - num=FLOOR(REAL(num2)/np) - IF (myid.LT.num2-(num2/np)*np) num=num+1 - k=1 - DO i=myid+1,num2,np - eig(k)=eig2(i) - k=k+1 - ENDDO - DEALLOCATE(eig2) - ! - ! Redistribute eigvec from ScaLAPACK distribution to each process - ! having all eigenvectors corresponding to his eigenvalues as above - ! - if (hamovlp%l_real) THEN - CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,zmat%z_r,eigvec_r) - ELSE - CALL subredist2(m,num2,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb,achi_c=zmat%z_c,asca_c=eigvec_c) - end if - ! - !DEALLOCATE ( eigvec) - DEALLOCATE ( iblacsnums ) - DEALLOCATE ( ihelp ) - DEALLOCATE ( iusermap ) - - END SUBROUTINE chani -end MODULE m_chani diff --git a/diagonalization/franza.F90 b/diagonalization/franza.F90 deleted file mode 100644 index 41044182e4118eda11f22790c78a9f5b82c6817f..0000000000000000000000000000000000000000 --- a/diagonalization/franza.F90 +++ /dev/null @@ -1,385 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -MODULE m_franza -CONTAINS - SUBROUTINE franza(nbasfcn,neigd, nsize,& - l_zref,l_J,matind,nred,gw,eig,ne,hamOvlp,zMat) - !*********************************************************************** - ! - ! solves the secular equation a*z=eig*b*z - ! - ! nbase,nval array dimensions as declared by calling routin - ! unchanged on exit - ! nsize actual dimension of a, b - ! a on entry: hamiltonian, lower triangle stored row-wise - ! = upper triangle stored column-wise - ! on exit : destroyed - ! b on entry: overlap matrix, lower tr. stored row-wise - ! on exit : destroyed(overwritten by cholesky factor) - ! eig on entry: ---- - ! on exit : eigenvalues in ascending rder - ! z on exit : corresponding eigenvectors - ! ne on exit : number of eigenvalues/vectors actually found - ! nblw on exit : number of eigenvalues below e(1) - ! work,iwork temporary storage of dimensions - ! nrtmp .ge. 8*nv - ! nitmp .ge. 5*nv - ! - ! This subroutine was written by Gustav Bihlmayer - ! IFF, July 1996 - ! - !*********************************************************************** - ! lapack routines inserted for the cray - ! please note, that if too many ev's are in your interval, you'll get - ! a memory fault. Then change dspevx('V','V' ... to dspevx('V','A' - ! to see what's going on .... - !*********************************************************************** - USE m_juDFT - USE m_types -#include"cpp_double.h" - IMPLICIT NONE - ! .. - ! .. Scalar Arguments .. - INTEGER, INTENT (IN) :: nbasfcn,neigd - INTEGER, INTENT (INOUT) :: nsize - INTEGER, INTENT (OUT) :: ne - INTEGER, INTENT (IN) :: nred,gw - LOGICAL, INTENT (IN) :: l_zref,l_J - ! .. - ! .. Array Arguments .. - INTEGER, INTENT (IN) :: matind(nbasfcn,2) - REAL, INTENT (OUT) :: eig(neigd) - TYPE(t_hamOvlp),INTENT(INOUT) :: hamOvlp - TYPE(t_zMat),INTENT(INOUT) :: zMat - - ! .. - ! .. Local Scalars .. - REAL toler,sq2i - INTEGER j,i1,i2,j1,j2,k,nsym,jsym,matsz - INTEGER info,i,iu,ne_a,nsize_a,n - LOGICAL addstate - ! .. - ! .. Local Arrays - REAL, ALLOCATABLE :: aa_r(:), bb_r(:) - REAL, ALLOCATABLE :: zz_r(:,:) - COMPLEX, ALLOCATABLE :: aa_c(:), bb_c(:), cwork(:) - COMPLEX, ALLOCATABLE :: zz_c(:,:) - - REAL, ALLOCATABLE :: work(:), etemp(:) - INTEGER, ALLOCATABLE :: iwork(:),ifail(:) - LOGICAL sort(nbasfcn) - REAL:: lb,ub - - !to select real/complex data - LOGICAL:: l_real - l_real=zMat%l_real -! IF (l_real.AND.PRESENT(a_c)) CALL juDFT_error("BUG in franza, call either with real OR complex data") - IF (l_real) THEN - ALLOCATE(zz_r(nbasfcn,neigd+1)) - ELSE - ALLOCATE(zz_c(nbasfcn,neigd+1)) - ENDIF - nsym=1 - IF (l_zref) THEN - ! - ! separate H and S matrix in symmetric (aa,bb) and antisymmetric (a,b) part - ! - IF (l_real) THEN - ALLOCATE (aa_r(nred*(nred+1)/2),bb_r(nred*(nred+1)/2)) - ELSE - ALLOCATE (aa_c(nred*(nred+1)/2),bb_c(nred*(nred+1)/2)) - END IF - - nsym=2 - i2=0 - j2=0 - DO i=1,nred - DO j=1,i - i1=(matind(i,1)-1)*matind(i,1)/2+matind(j,1) - j1=(matind(i,1)-1)*matind(i,1)/2+matind(j,2) - i2=i2+1 - IF (l_real) THEN - aa_r(i2)=hamOvlp%a_r(i1)+hamOvlp%a_r(j1) - bb_r(i2)=hamOvlp%b_r(i1)+hamOvlp%b_r(j1) - ELSE - aa_c(i2)=hamOvlp%a_c(i1)+hamOvlp%a_c(j1) - bb_c(i2)=hamOvlp%b_c(i1)+hamOvlp%b_c(j1) - END IF - IF ((matind(i,1).NE.matind(i,2)).AND.(matind(j,1).NE.matind(j,2))) THEN - j2=j2+1 - IF (l_real) THEN - hamOvlp%a_r(j2)=hamOvlp%a_r(i1)-hamOvlp%a_r(j1) - hamOvlp%b_r(j2)=hamOvlp%b_r(i1)-hamOvlp%b_r(j1) - ELSE - hamOvlp%a_c(j2)=hamOvlp%a_c(i1)-hamOvlp%a_c(j1) - hamOvlp%b_c(j2)=hamOvlp%b_c(i1)-hamOvlp%b_c(j1) - END IF - ENDIF - ENDDO - ENDDO - nsize = nsize - nred - ENDIF - DO jsym=1,nsym - - IF (jsym.EQ.2) THEN - ! - ! second time calculate symmetric part and store the antisym. EV's and EW's - ! - matsz=(nred+1)*nred/2 - IF (l_real) THEN - CALL CPP_BLAS_scopy(matsz,aa_r,1,hamOvlp%a_r,1) - CALL CPP_BLAS_scopy(matsz,bb_r,1,hamOvlp%b_r,1) - ELSE - CALL CPP_BLAS_ccopy(matsz,aa_c,1,hamOvlp%a_c,1) - CALL CPP_BLAS_ccopy(matsz,bb_c,1,hamOvlp%b_c,1) - ENDIF - - ne_a=ne - k=1 - nsize_a=nsize - DO i=1,ne - IF (l_real) THEN - aa_r(i)=eig(i) - ELSE - aa_c(i)=CMPLX(eig(i),0.0) - ENDIF - IF (l_real) THEN - DO j=1,nsize - bb_r(k)=zMat%z_r(j,i) - k=k+1 - ENDDO - ELSE - DO j=1,nsize - bb_c(k)=zMat%z_c(j,i) - k=k+1 - ENDDO - ENDIF - ENDDO - nsize=nred - ENDIF - !-gu - ! --> start with Cholesky factorization of b ( so that b = l * l^t) - ! --> b is overwritten by l - ! - IF (l_real) THEN - CALL CPP_LAPACK_spptrf('U',nsize,hamOvlp%b_r,info) - ELSE - CALL CPP_LAPACK_cpptrf('U',nsize,hamOvlp%b_c,info) - ENDIF - ! - IF (info.NE.0) THEN - WRITE (*,*) 'Error in cpptrf/spptrf: info =',info - CALL juDFT_error("Diagonalization failed",calledby="franza") - ENDIF - ! - ! --> now reduce a * z = eig * b * z to the standard form a' * z' = eig * z' - ! --> where a' = (l)^-1 * a * (l^t)^-1 and z' = l^t * z - ! - IF (l_real) THEN - CALL CPP_LAPACK_sspgst(1,'U',nsize,hamOvlp%a_r,hamOvlp%b_r,info) - ELSE - CALL CPP_LAPACK_chpgst(1,'U',nsize,hamOvlp%a_c,hamOvlp%b_c,info) - ENDIF - ! - IF (info.NE.0) THEN - WRITE (6,*) 'Error in chpgst/sspgst: info =',info - CALL juDFT_error("Diagonalization failed",calledby="franza") - ENDIF - ! - ! --> solve a' * z' = eig * z' for eigenvalues eig between lb und ub - ! - iu = MIN(neigd,nsize) - IF (l_zref) iu = (neigd-1)/2 - ! toler=2.0*CPP_LAPACK_slamch('S') - toler=2.0*TINY(sq2i) - ALLOCATE ( work(8*nbasfcn),iwork(5*nbasfcn),ifail(nbasfcn) ) - ALLOCATE ( etemp(nbasfcn) ) - addstate = gw.NE.0.AND.iu.LT.nsize ! add one state, - IF(addstate) iu = iu + 1 ! see below (CF) - - IF (l_real) THEN - zz_r=0.0 - IF(l_J)THEN - CALL CPP_LAPACK_sspevx('N','I','U',nsize,hamOvlp%a_r,lb,ub,1,iu,toler,ne, etemp,zz_r,nbasfcn,work,iwork,ifail,info) - ELSE - CALL CPP_LAPACK_sspevx('V','I','U',nsize,hamOvlp%a_r,lb,ub,1,iu,toler,ne, etemp,zz_r,nbasfcn,work,iwork,ifail,info) - ENDIF - ELSE - zz_c=0.0 - ALLOCATE ( cwork(2*nbasfcn) ) - IF(l_J)THEN - CALL CPP_LAPACK_chpevx('N','I','U',nsize,hamOvlp%a_c,lb,ub,1,iu,toler,ne, etemp,zz_c,nbasfcn,cwork,work,iwork,ifail,info) - ELSE - CALL CPP_LAPACK_chpevx('V','I','U',nsize,hamOvlp%a_c,lb,ub,1,iu,toler,ne, etemp,zz_c,nbasfcn,cwork,work,iwork,ifail,info) - ENDIF - DEALLOCATE ( cwork ) - ENDIF - IF((info.EQ.0).AND.(ne.NE.iu)) THEN - WRITE(*,*) '' - WRITE(*,*) 'Call to LAPACK routine sspevx/dspevx/chpevx/zhpevx resulted in problem:' - WRITE(*,*) 'The library does not report an error but the number of found eigenvalues' - WRITE(*,*) 'does not agree with the target number of eigenvalues to be calculated.' - WRITE(*,*) 'This is a bug in either fleur or the linked LAPACK library.' - WRITE(*,*) '' - CALL juDFT_error("sspevx/dspevx/chpevx/zhpevx failed without explicit errorcode.",calledby="franza") - END IF - IF(addstate) THEN ! cut topmost subspace of degenerate states to avoid symmetry breaking (CF) - iu = ne - ne = ne - 1 - DO WHILE(ABS(etemp(ne)-etemp(iu)).LT.1d-8) - ne = ne - 1 - ENDDO - ENDIF - IF (l_real) THEN - zMat%z_r(:,:ne) = zz_r(:,:ne) - ELSE - zMat%z_c(:,:ne) = zz_c(:,:ne) - END IF - ! - IF (ne.GT.neigd) THEN - WRITE(6,*) 'ne=',ne,' > neigd' - CALL juDFT_error("ne.GT.neigd",calledby="franza") - ENDIF - eig(:ne) = etemp(:ne) - ! - IF(l_J) THEN - DEALLOCATE ( work,iwork,ifail,etemp ) - ELSE - IF (info.NE.0) THEN - WRITE (6,*) 'Error in chpexvx/sspevx: info =',info - WRITE (6,*) 'The following eigenvectors did not converge:' - WRITE (6,'(30i5)') (ifail(i),i=1,ne) - CALL juDFT_error("Diagonalization failed",calledby="franza") - ENDIF - DEALLOCATE ( work,iwork,ifail,etemp ) - ! - ! --> recover the generalized eigenvectors z by solving z' = l^t * z - ! - IF (l_real) THEN - CALL CPP_LAPACK_stptrs('U','N','N',nsize,ne,hamOvlp%b_r,zMat%z_r,nbasfcn,info) - ELSE - CALL CPP_LAPACK_ctptrs('U','N','N',nsize,ne,hamOvlp%b_c,zMat%z_c,nbasfcn,info) - ENDIF - ! - IF (info.NE.0) THEN - WRITE (6,*) 'Error in c/stptrs: info =',info - CALL juDFT_error("Diagonalization failed",calledby="franza") - ENDIF - ENDIF !l_J - !+gu - ENDDO - ! - ! now collect symmetric and antisym. EW's and EV's and sort - ! - IF (l_zref) THEN - k=1 - DO i=1,ne - IF (l_real) THEN - DO j=1,nsize - hamOvlp%b_r(k)=zMat%z_r(j,i) - k=k+1 - ENDDO - ELSE - DO j=1,nsize - hamOvlp%b_c(k)=zMat%z_c(j,i) - k=k+1 - ENDDO - END IF - ENDDO - ! - ! prepare sort-array: even=.true., odd=.false. - ! - i=1 - j=1 - eig(ne+1) = 99.9e9 - IF (l_real) THEN - aa_r(ne_a+1) = 99.9e9 - ELSE - aa_c(ne_a+1) = CMPLX(99.9e9,0.0) - ENDIF - DO k=1,ne+ne_a - IF (l_real) THEN - IF (eig(i).LT.REAL(aa_r(j))) THEN - sort(k)=.TRUE. - i=i+1 - ELSE - sort(k)=.FALSE. - j=j+1 - ENDIF - ELSE - IF (eig(i).LT.REAL(aa_c(j))) THEN - sort(k)=.TRUE. - i=i+1 - ELSE - sort(k)=.FALSE. - j=j+1 - ENDIF - ENDIF - ENDDO - ! - ! sort EW's and EV's - ! - i=ne - j=ne_a - nsize = nsize + nsize_a - ne = ne + ne_a - sq2i=1.0/SQRT(2.0) - DO k=ne,1,-1 - DO n=1,nsize - IF (l_real) THEN - zMat%z_r(n,k)=0.0 - ELSE - zMat%z_c(n,k)=CMPLX(0.0,0.0) - ENDIF - ENDDO - IF (sort(k)) THEN - eig(k)=eig(i) - i1=nred * (i-1) - DO n=1,nred - i1=i1+1 - IF (l_real) THEN - zMat%z_r(matind(n,1),k) = zMat%z_r(matind(n,1),k)+hamOvlp%b_r(i1)*sq2i - zMat%z_r(matind(n,2),k) = zMat%z_r(matind(n,2),k)+hamOvlp%b_r(i1)*sq2i - ELSE - zMat%z_c(matind(n,1),k) = zMat%z_c(matind(n,1),k)+hamOvlp%b_c(i1)*CMPLX(sq2i,0.0) - zMat%z_c(matind(n,2),k) = zMat%z_c(matind(n,2),k)+hamOvlp%b_c(i1)*CMPLX(sq2i,0.0) - ENDIF - ENDDO - i=i-1 - ELSE - IF (l_real) THEN - eig(k)=aa_r(j) - ELSE - eig(k)=REAL(aa_c(j)) - ENDIF - j1=nsize_a * (j-1) - DO n=1,nred - IF (matind(n,1).NE.matind(n,2)) THEN - j1=j1+1 - IF (l_real) THEN - zMat%z_r(matind(n,1),k) = zMat%z_r(matind(n,1),k)+bb_r(j1)*sq2i - zMat%z_r(matind(n,2),k) = zMat%z_r(matind(n,2),k)-bb_r(j1)*sq2i - ELSE - zMat%z_c(matind(n,1),k) = zMat%z_c(matind(n,1),k)+bb_c(j1) *CMPLX(sq2i,0.0) - zMat%z_c(matind(n,2),k) = zMat%z_c(matind(n,2),k)-bb_c(j1) *CMPLX(sq2i,0.0) - ENDIF - ELSE - IF (l_real) THEN - zMat%z_r(matind(n,1),k) = 0.0 - ELSE - zMat%z_c(matind(n,1),k) = CMPLX(0.0,0.0) - ENDIF - ENDIF - ENDDO - j=j-1 - ENDIF - ENDDO - - ENDIF - !-gu - END SUBROUTINE franza -END diff --git a/diagonalization/geneigbrobl.F90 b/diagonalization/geneigbrobl.F90 deleted file mode 100644 index 43643d4e9673f9bfb68d9a4e9c2e09b48d05661d..0000000000000000000000000000000000000000 --- a/diagonalization/geneigbrobl.F90 +++ /dev/null @@ -1,173 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -MODULE m_geneigprobl - USE m_juDFT - !********************************************************** - ! Solve the generalized eigenvalue problem - ! Frank Freimuth, November 2006 - !********************************************************** -CONTAINS - SUBROUTINE geneigprobl(nbasfcn, nsize,neigd,l_J,eig,ne,a_r,b_r,z_r,a_c,b_c,z_c) -#include"cpp_double.h" - USE m_packed_to_full - IMPLICIT NONE - - ! ... Arguments ... - - INTEGER, INTENT (IN) :: nbasfcn - INTEGER, INTENT (IN) :: neigd - INTEGER, INTENT (IN) :: nsize - LOGICAL, INTENT (IN) :: l_J - - REAL, INTENT(OUT) :: eig(:) - INTEGER, INTENT(OUT) :: ne - - REAL,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: a_r(:),b_r(:) - REAL,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: z_r(:,:) - COMPLEX,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: a_c(:),b_c(:) - COMPLEX,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: z_c(:,:) - - - ! ... Local Variables .. - - INTEGER iind,ind1,ind2,info,lwork,liwork,lrwork,err,i - INTEGER sizez,iu - REAL :: lb,ub - ! 'sizez' is needed, as some compilers sometimes produce errors, - ! if the size command is used directly as a lapack argument. - - REAL toler, eigTemp(nsize) - - REAL, ALLOCATABLE :: work(:) - INTEGER, ALLOCATABLE :: iwork(:),isuppz(:) - REAL, ALLOCATABLE :: largea_r(:,:),largeb_r(:,:) - COMPLEX, ALLOCATABLE :: largea_c(:,:),largeb_c(:,:) - COMPLEX,ALLOCATABLE :: cwork(:) - - LOGICAL :: l_real - - l_real=PRESENT(a_r) - - - !********************************** - !expand from packed to full storage: full storage lapack-routines - !are faster than the packed lapack-routines. - !********************************** - !hamiltonian - IF (l_real) THEN - call packed_to_full(nsize,a_r,largea_r) - DEALLOCATE (a_r) - call packed_to_full(nsize,b_r,largeb_r) - DEALLOCATE (b_r) - ELSE - call packed_to_full(nsize,a_c,largea_c) - DEALLOCATE (a_c) - call packed_to_full(nsize,b_c,largeb_c) - DEALLOCATE (b_c) - ENDIF - - - - IF (l_real) then - CALL CPP_LAPACK_spotrf('U',nsize,largeb_r,nsize,info) - IF (info /= 0) CALL juDFT_error("error in spotrf",calledby ="geneigprobl") - - CALL CPP_LAPACK_ssygst(1,'U',nsize,largea_r,nsize,largeb_r,nsize,info) - IF (info /= 0) CALL juDFT_error("error in ssygst",calledby ="geneigprobl") - - toler = 2.0*TINY(toler) - liwork = 10*nsize - ALLOCATE ( iwork(liwork), stat=err ) - IF (err/=0) CALL juDFT_error("error allocating iwork",calledby ="geneigprobl") - - lwork = 26*nsize - ALLOCATE ( work(lwork), stat=err ) - IF (err/=0) CALL juDFT_error(" error allocating work",calledby ="geneigprobl") - ALLOCATE ( isuppz(2*nsize), stat=err ) - IF (err /= 0) CALL juDFT_error("error allocating isuppz",calledby ="geneigprobl") - IF (ALLOCATED(z_r)) THEN - IF (.NOT.(SIZE(z_r,1)==nbasfcn.AND.SIZE(z_r,2)==neigd)) DEALLOCATE(z_r) - ENDIF - IF (.NOT.ALLOCATED(z_r)) THEN - ALLOCATE ( z_r(nbasfcn,neigd), stat=err ) - IF (err/=0) THEN - WRITE(*,*) nbasfcn,neigd,err - CALL juDFT_error("error allocating z",calledby ="geneigprobl") - ENDIF - ENDIF - sizez= SIZE(z_r,1) - iu = MIN(nsize,neigd) - IF (l_J) THEN - CALL CPP_LAPACK_ssyevr('N','I','U',nsize,largea_r, nsize,lb,ub,1,iu,toler,ne,eigTemp,z_r,& - sizez,isuppz,work,lwork,iwork,liwork,info) - ELSE - CALL CPP_LAPACK_ssyevr('V','I','U',nsize,largea_r,nsize,lb,ub,1,iu,toler,ne,eigTemp,z_r,& - sizez,isuppz,work,lwork,iwork,liwork,info) - ENDIF - IF (info /= 0) CALL juDFT_error("error in ssyevr",calledby ="geneigprobl") - DEALLOCATE (isuppz,work,iwork) - - CALL CPP_LAPACK_strtrs('U','N','N',nsize,ne,largeb_r, nsize,z_r,sizez,info) - IF (info /= 0) CALL juDFT_error("error in strtrs",calledby ="geneigprobl") - ELSE - - CALL CPP_LAPACK_cpotrf('U',nsize,largeb_c,nsize,info) - IF (info /= 0) CALL juDFT_error("error in cpotrf",calledby ="geneigprobl") - - CALL CPP_LAPACK_chegst(1,'U',nsize,largea_c,nsize,largeb_c,nsize,info) - IF (info /= 0) CALL juDFT_error(" error in chegst",calledby ="geneigprobl") - - toler = 2.0*TINY(toler) - liwork = 50*nsize - ALLOCATE ( iwork(liwork), stat=err ) - IF (err/=0) CALL juDFT_error("error allocating iwork",calledby ="geneigprobl") - - lwork = 20*nsize - ALLOCATE( cwork(lwork), stat=err ) - IF (err/=0) CALL juDFT_error("error allocating cwork",calledby ="geneigprobl") - ALLOCATE( isuppz(10*nsize), stat=err ) - IF (err/=0) CALL juDFT_error("error allocating isuppz",calledby ="geneigprobl") - - lrwork = 84*nsize - ALLOCATE (work(lrwork), stat=err ) - IF (err/=0) CALL juDFT_error(" error allocating work",calledby ="geneigprobl") - IF (ALLOCATED(z_c)) THEN - IF (.NOT.(SIZE(z_c,1)==nbasfcn.AND.SIZE(z_c,2)==neigd)) DEALLOCATE(z_c) - ENDIF - IF (.NOT.ALLOCATED(z_c)) THEN - ALLOCATE ( z_c(nbasfcn,neigd), stat=err ) - IF (err/=0) THEN - WRITE(*,*) nbasfcn,neigd,err - CALL juDFT_error("error allocating z",calledby ="geneigprobl") - ENDIF - ENDIF - sizez= SIZE(z_c,1) - iu = MIN(nsize,neigd) - IF (l_J) THEN - CALL CPP_LAPACK_cheevr('N','I','U',nsize,largea_c, nsize,lb,ub,1,iu,toler,ne,eigTemp,z_c,& - sizez,isuppz,cwork,lwork,work,lrwork,iwork,liwork,info) - ELSE - CALL CPP_LAPACK_cheevr('V','I','U',nsize,largea_c, nsize,lb,ub,1,iu,toler,ne,eigTemp,z_c,& - sizez,isuppz,cwork,lwork,work,lrwork,iwork,liwork,info) - ENDIF - IF (info /= 0) CALL juDFT_error("error in cheevr",calledby ="geneigprobl") - DEALLOCATE ( isuppz ) - DEALLOCATE ( work ) - DEALLOCATE ( iwork ) - DEALLOCATE ( cwork ) - - CALL CPP_LAPACK_ctrtrs('U','N','N',nsize,ne,largeb_c, nsize,z_c,sizez,info) - IF (info /= 0) CALL juDFT_error("error in ctrtrs",calledby ="geneigprobl") - ENDIF - - DO i = 1, neigd - eig(i) = eigTemp(i) - END DO - -END SUBROUTINE geneigprobl -END MODULE m_geneigprobl - diff --git a/diagonalization/locrectify.F90 b/diagonalization/locrectify.F90 deleted file mode 100644 index b86ef68c4bb4f70319e0cc3752d71e0fa9fc08b9..0000000000000000000000000000000000000000 --- a/diagonalization/locrectify.F90 +++ /dev/null @@ -1,445 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -!************************************************************* -! Find linear combinations of local orbitals that are -! eigenfunctions of the z-reflection operation. -! Used to exploit z-reflection symmetry when solving -! the eigenvalue problem (film calculations). -! Frank Freimuth, January 2007 -!************************************************************* -#include"cpp_double.h" -MODULE m_locrectify - USE m_juDFT - CONTAINS - SUBROUTINE locrectify(jsp,input,lapw,bkpt,atoms, kveclo, sym,cell,& - kindlocrec,evenlocs,oddlocs,evenlocindex, oddlocindex,locrec_r,locrec_c) - - USE m_ylm - USE m_constants,ONLY:tpi_const - USE m_loccoeff - USE m_types - IMPLICIT NONE - INTEGER,INTENT(IN) :: jsp - TYPE(t_input),INTENT(IN) :: input - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_lapw),INTENT(IN) :: lapw - REAL,INTENT(in)::bkpt(3) - INTEGER,INTENT(in)::kveclo(atoms%nlotot) - - REAL,OPTIONAL,INTENT(out)::locrec_r(atoms%nlotot,atoms%nlotot) - COMPLEX,OPTIONAL,INTENT(out)::locrec_c(atoms%nlotot,atoms%nlotot) - - LOGICAL,INTENT(out)::kindlocrec(atoms%nlotot) - INTEGER,INTENT(out)::evenlocs - INTEGER,INTENT(out)::oddlocs - INTEGER,INTENT(out)::evenlocindex(atoms%nlotot) - INTEGER,INTENT(out)::oddlocindex(atoms%nlotot) - - INTEGER numorec - INTEGER invatom,zrefatom - LOGICAL l_planetwo,l_spacetwo,l_spacefour,l_spacetwoinv - LOGICAL l_crosstwo - INTEGER pos,loc,nn,angmom,bas,locm,loopi,j,i - INTEGER lm,locmm,lmm - INTEGER ind,vec,natom,nap,inap,locnum,ll1 - INTEGER ipiv(28),info - COMPLEX loccoffs(28,28),reccoffs(28) - COMPLEX ylm(16),phase - REAL fk(3),tmk,fkp(3) - ! complex coeffi(1:28,1:28,0:13) transferred to mod_loccoeff.f - ! logical coffkind(1:28,0:13) import with USE m_loccoeff - INTEGER jump(atoms%nat) - INTEGER mqnum,nalo,naup,basindex(atoms%nat,atoms%nlod) - INTEGER zrefnap,zrefinap,coffindex - LOGICAL l_real,l_invsatswap - - ! ci=cmplx(0.0,1.0) transferred to mod_loccoeff.f - - l_real=PRESENT(locrec_r) - - jump(1:atoms%nat)=0 - IF (l_real) THEN - locrec_r=0.0 - ELSE - locrec_c=0.0 - ENDIF - numorec=0 - - !*********************************************** - ! basindex holds first index of subset of locs - !*********************************************** - bas=1 - nalo=1 - DO ind=1,atoms%ntype - naup=nalo+atoms%neq(ind)-1 - DO natom=nalo,naup - IF(atoms%invsat(natom).EQ.2)CYCLE - DO loc=1,atoms%nlo(ind) - basindex(natom,loc)=bas - locnum=(atoms%invsat(natom)+1)*(2*atoms%llo(loc,ind)+1) - bas=bas+locnum - ENDDO !loc - ENDDO !natom - nalo=naup+1 - ENDDO !ind - - bas=0 - nalo=1 - DO ind=1,atoms%ntype - naup=nalo+atoms%neq(ind)-1 - DO natom=nalo,naup - IF(atoms%invsat(natom).EQ.2)CYCLE - IF(jump(natom).NE.0)THEN - bas=bas+jump(natom) - CYCLE - ENDIF - - !******************************************************* - ! find out by what kind of symmetry the atoms - ! are interrelated - !******************************************************* - !two atoms related by inversion symmetry are in xy-plane - l_planetwo=.FALSE. - !two atoms are related by z-reflection, but not by inversion - l_spacetwo=.FALSE. - !two atoms are related by z-reflection and by inversion - l_spacetwoinv=.FALSE. - !four atoms are entangled, by both inversion and z-reflection - l_spacefour=.FALSE. - !two atoms are related by inversion, but not by z-reflection - l_crosstwo=.FALSE. - !order of locs corresponds to setup of coeffi-array => l_invsatswap=.false. - !order of locs is reversed with respect to coeffi-array => l_invsatswap=.true. - l_invsatswap=.FALSE. - IF((ABS(ABS(atoms%taual(3,natom))-0.5).LT.1e-6).AND..NOT.input%film)THEN !atom on face - IF(atoms%invsat(natom).EQ.1)THEN ! two atoms - invatom=sym%invsatnr(natom) - l_planetwo=.TRUE. - ENDIF - ELSEIF(ABS(atoms%taual(3,natom)).GT.1e-6) THEN !atom not in xy-plane - DO zrefatom=nalo,naup !find z-reflection image - IF(ALL(ABS(atoms%taual(1:2,natom)-atoms%taual(1:2,zrefatom))-& - NINT(ABS(atoms%taual(1:2,natom)-atoms%taual(1:2,zrefatom))).LT.1e-6)) THEN - IF(ABS(atoms%taual(3,natom)+atoms%taual(3,zrefatom)).LT.1e-6) GOTO 543 - ENDIF - ENDDO !zrefatom - if (l_real) THEN - !if there is no z-reflected counterpart there must be an inversion - !counterpart instead - IF(atoms%invsat(natom).EQ.1)THEN - zrefatom=0 !this switches on l_crosstwo later on - ELSE - CALL juDFT_error("nozref-invsat",calledby ="locrectify") - ENDIF - else - !in the complex version there must be a z-reflected counterpart - CALL juDFT_error("zrefatom not found",calledby ="locrectify") - endif -543 CONTINUE - - IF(atoms%invsat(natom).EQ.0) THEN !no inversion-image - l_spacetwo=.TRUE. - if (l_real) CALL juDFT_error("spacetwo & INVERSION",calledby ="locrectify") - ELSE !inversion-image - invatom=sym%invsatnr(natom) - IF(zrefatom.EQ.0)THEN - l_crosstwo=.TRUE. - ELSEIF(invatom.EQ.zrefatom)THEN !inversion = z-reflection - l_spacetwoinv=.TRUE. - ELSE !inversion /= z-reflection - if (l_real) THEN - l_spacefour=.TRUE. !four entangled locs - else - CALL juDFT_error("spacefour & no INVERSION",calledby ="locrectify") - endif - IF(atoms%invsat(zrefatom).EQ.1)THEN - - ELSEIF(atoms%invsat(zrefatom).EQ.2)THEN - l_invsatswap=.TRUE. - DO info=1,atoms%nat - IF(sym%invsatnr(info).EQ.zrefatom)GOTO 723 - ENDDO - CALL juDFT_error("no atoms for zref",calledby ="locrectify") -723 CONTINUE - zrefatom=info - IF(atoms%invsat(zrefatom)/=1) CALL juDFT_error("atoms%invsat" ,calledby ="locrectify") - IF(.NOT.ALL(ABS(atoms%taual(1:2,invatom)- atoms%taual(1:2,zrefatom))<1e-6)) & - CALL juDFT_error ("invsat-zref2",calledby ="locrectify") - IF(.NOT.ABS(atoms%taual(3,invatom)+atoms%taual(3,zrefatom)) <1e-6) & - CALL juDFT_error("invsat-zref3", calledby="locrectify") - - ELSEIF(atoms%invsat(zrefatom).EQ.0)THEN - CALL juDFT_error("spacefour: atoms-zref",calledby ="locrectify") - ENDIF - ENDIF - ENDIF - ELSE !atom lies in the xy-plane - IF(atoms%invsat(natom).EQ.1)THEN ! two atoms - invatom=sym%invsatnr(natom) - l_planetwo=.TRUE. - ENDIF - ENDIF !symmetry partners - - DO loc=1,atoms%nlo(ind) - angmom=atoms%llo(loc,ind) - IF(angmom.EQ.0)THEN - IF(l_planetwo)THEN - kindlocrec(numorec+1)=.TRUE. - kindlocrec(numorec+2)=.TRUE. - IF (l_real) THEN - locrec_r(bas+1,numorec+1)=1.0 - locrec_r(bas+2,numorec+2)=1.0 - ELSE - locrec_c(bas+1,numorec+1)=1.0 - locrec_c(bas+2,numorec+2)=1.0 - ENDIF - bas=bas+2 - numorec=numorec+2 - CYCLE - ELSEIF(l_spacetwoinv.OR.l_crosstwo)THEN - coffindex=0 - locnum=2 - ELSEIF(l_spacetwo) THEN - locnum=1 - coffindex=0 - ELSEIF(l_spacefour)THEN - locnum=2 - coffindex=6 - ELSE - kindlocrec(numorec+1)=.TRUE. - IF (l_real) THEN - locrec_r(bas+1,numorec+1)=1.0 - ELSE - locrec_c(bas+1,numorec+1)=1.0 - END IF - bas=bas+1 - numorec=numorec+1 - CYCLE - ENDIF !l_planetwo,l_crosstwo,.... - ELSE - IF(l_planetwo)THEN - locnum=2*(2*angmom+1) - IF(angmom.EQ.1)THEN - coffindex=3 - ELSEIF(angmom.EQ.2)THEN - coffindex=9 - ELSEIF(angmom.EQ.3)THEN - coffindex=12 - ELSE - CALL juDFT_error("angmom1",calledby ="locrectify") - ENDIF - ELSEIF(l_spacetwoinv.OR.l_crosstwo)THEN - locnum=2*(2*angmom+1) - IF(angmom.EQ.1)THEN - coffindex=4 - ELSEIF(angmom.EQ.2)THEN - coffindex=8 - ELSEIF(angmom.EQ.3)THEN - coffindex=11 - ELSE - CALL juDFT_error("angmom5",calledby ="locrectify") - ENDIF - ELSEIF(l_spacetwo) THEN - locnum=2*angmom+1 - IF(angmom.EQ.1)THEN - coffindex=4 - ELSEIF(angmom.EQ.2)THEN - coffindex=8 - ELSEIF(angmom.EQ.3)THEN - coffindex=11 - ELSE - CALL juDFT_error("angmom2",calledby="locrectify") - ENDIF - ELSEIF(l_spacefour) THEN - locnum=2*(2*angmom+1) - IF(angmom.EQ.1)THEN - coffindex=7 - ELSEIF(angmom.EQ.2)THEN - coffindex=10 - ELSEIF(angmom.EQ.3)THEN - coffindex=13 - ELSE - CALL juDFT_error("angmom3",calledby ="locrectify") - ENDIF - ELSE - locnum=2*angmom+1 - IF(angmom.EQ.1)THEN - coffindex=1 - ELSEIF(angmom.EQ.2)THEN - coffindex=2 - ELSEIF(angmom.EQ.3)THEN - coffindex=5 - ELSE - CALL juDFT_error("angmom",calledby ="locrectify") - ENDIF - ENDIF! l_planetwo - ENDIF !angmom - - mqnum=2*angmom+1 - ll1=angmom*angmom-1 - loccoffs(:,:)=CMPLX(0.0,0.0) - !************************************************** - ! Write the expansion of the locs - ! in terms of spherical harmonics into the - ! loccoffs-array. - ! First index of loccoffs: m-quantum-number. - ! Second index of loccoffs: index of loc. - !************************************************** - DO locm=1,locnum - pos=bas+locm - vec=kveclo(pos) - !The vector of the planewave onto which the loc is matched. - fk(:)=bkpt(:)+(/lapw%k1(vec,jsp),lapw%k2(vec,jsp),lapw%k3(vec,jsp)/) - - tmk=tpi_const*DOT_PRODUCT(lapw%k1(:,jsp),atoms%taual(:,natom)) - phase = CMPLX(COS(tmk),SIN(tmk)) - fkp=MATMUL(fk,cell%bmat) - CALL ylm4(3,fkp,ylm) - DO locmm=1,mqnum - lmm=ll1+locmm - loccoffs(locmm,locm)=ci**angmom*phase*CONJG(ylm(lmm+1)) - ENDDO - ENDDO !locm - !************************************************************* - ! If locs are entangled by symmetry, the number of locs - ! that have to be treated at the same time is larger. - ! => l_spacetwo: two locs, which are NOT in the - ! xy-plane and which are not related by inversion, are - ! entangled by z-reflexion. - ! => l_spacefour: two pairs the partners of which are - ! connected by inversion are interrelated by z-reflection. - !************************************************************* - IF(l_spacetwo.OR.l_spacefour)THEN - DO locm=1,locnum - pos=basindex(zrefatom,loc)-1+locm - vec=kveclo(pos) - fk(:)=bkpt(:)+(/lapw%k1(vec,jsp),lapw%k2(vec,jsp),lapw%k3(vec,jsp)/) - - tmk=tpi_const*DOT_PRODUCT(lapw%k1(:,jsp),atoms%taual(:,zrefatom)) - phase = CMPLX(COS(tmk),SIN(tmk)) - fkp=MATMUL(fk,cell%bmat) - CALL ylm4(3,fkp,ylm) - DO locmm=1,mqnum - lmm=ll1+locmm - loccoffs(locnum+locmm,locnum+locm)= ci**angmom*phase*CONJG(ylm(lmm+1)) - ENDDO - ENDDO !locm - jump(zrefatom)=jump(zrefatom)+locnum - IF(l_spacetwo)locnum=locnum*2 - ENDIF !l_spacetwo - !************************************************************* - ! If locs are entangled by symmetry, loccoffs has a larger - ! number of components for a given loc. - ! l_planetwo => two locs in xy-plane are entangled - ! by inversion symmetry. - ! l_spacetwoinv => two z-reflexion-symmetric atoms - ! are entangled also by inversion. - ! l_crosstwo => two inversion-symmetric atoms - ! are not directly connected by z-reflexion. - ! l_spacefour => two pairs the partners of which are - ! connected by inversion are interrelated by z-reflection. - !************************************************************* - IF(l_planetwo.OR.l_crosstwo.OR.l_spacetwoinv .OR.l_spacefour) THEN - DO locm=1,locnum - DO locmm=1,mqnum - loccoffs(mqnum+locmm,locm)=(-1)**(locmm-1)* CONJG(loccoffs(mqnum+1-locmm,locm)) - ENDDO! locmm - ENDDO ! locm - ENDIF ! l_planetwo - - !****************************************************************** - ! If four atoms are entangled, something remains to be added - ! here: - !****************************************************************** - IF(l_spacefour)THEN - DO locm=1,locnum - DO locmm=1,mqnum - loccoffs(locnum+mqnum+locmm,locnum+locm)=(-1)**(locmm-1)*& - CONJG(loccoffs(locnum+mqnum+1-locmm,locnum+locm)) - ENDDO! locmm - ENDDO ! locm - locnum=locnum*2 - ENDIF - !****************************************************************** - ! Find linear combinations of locs that are eigenfunctions - ! of the z-reflection operation by solving linear system - ! of equations. Write the transformation into locrec-matrix. - !****************************************************************** - CALL CPP_LAPACK_cgetrf(locnum,locnum,loccoffs,28,ipiv,info) - IF(info /= 0) CALL juDFT_error("cgetrf",calledby ="locrectify") - DO loopi=1,locnum - numorec=numorec+1 - kindlocrec(numorec)=coffkind(loopi,coffindex) - IF(l_invsatswap)THEN - reccoffs(1:locnum/2)=coeffi(1:locnum/2,loopi,coffindex) - reccoffs(locnum/2+1:locnum/4*3)= coeffi(locnum/4*3+1:locnum,loopi,coffindex) - reccoffs(locnum/4*3+1:locnum)= coeffi(locnum/2+1:locnum/4*3,loopi,coffindex) - ELSE - reccoffs(1:locnum)=coeffi(1:locnum,loopi,coffindex) - ENDIF - CALL CPP_LAPACK_cgetrs('N',locnum,1,loccoffs,28,ipiv, reccoffs,28,info) - IF(info /= 0) CALL juDFT_error("cgetrs", calledby="locrectify") - IF (l_real) THEN - IF(ANY(ABS(AIMAG(reccoffs(1:locnum))).GT.1e-6)) THEN - !In the real version of the code, the hamiltonian and overlap - !matrices are real. Hence the transformation matrices that make - !the locs become eigenfunctions of the z-reflection operation - !have to be real also. Check this here. - PRINT*,"Sorry!" - PRINT*,"The transformation is not purely real." - PRINT*,"=> Something is wrong here." - CALL juDFT_error("unfortunately complex",calledby ="locrectify") - ENDIF - IF(l_spacefour)THEN - locrec_r(bas+1:bas+locnum/2,numorec)=REAL(reccoffs(:locnum/2)) - DO locmm=1,locnum/2 - locrec_r(basindex(zrefatom,loc)-1+locmm,numorec)= REAL(reccoffs(locmm+locnum/2)) - ENDDO - - ELSE - locrec_r(bas+1:bas+locnum,numorec)=REAL(reccoffs(:locnum)) - - ENDIF - ELSE - IF(l_spacetwo)THEN - locrec_c(bas+1:bas+locnum/2,numorec)=reccoffs(1:locnum/2) - locrec_c(basindex(zrefatom,loc):basindex(zrefatom,loc)+locnum/2-1,numorec)=& - reccoffs(1+locnum/2:locnum) - ELSE - locrec_c(bas+1:bas+locnum,numorec)=reccoffs(1:locnum) - ENDIF - ENDIF - ENDDO !loopi - IF(l_spacetwo) locnum=locnum/2 - IF(l_spacefour) locnum=locnum/2 - bas=bas+locnum - ENDDO !loc - ENDDO !natom - nalo=naup+1 - ENDDO !ind - - - !********************************************************* - ! Determine number of even and odd locs. - ! Prepare sort-arrays. - !********************************************************* - evenlocs=0 - oddlocs=0 - DO loopi=1,atoms%nlotot - IF(kindlocrec(loopi)) THEN - evenlocs=evenlocs+1 - evenlocindex(evenlocs)=loopi - ELSE - oddlocs=oddlocs+1 - oddlocindex(oddlocs)=loopi - ENDIF - ENDDO - END SUBROUTINE locrectify -END MODULE m_locrectify - diff --git a/diagonalization/mod_loccoeff.f90 b/diagonalization/mod_loccoeff.f90 deleted file mode 100644 index f8994d1b4ffec2c6a1ec26dc4e5fda7e965ffddf..0000000000000000000000000000000000000000 --- a/diagonalization/mod_loccoeff.f90 +++ /dev/null @@ -1,868 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -MODULE m_loccoeff - - ! coefficients used in locrectify: - - COMPLEX coeffi(1:28,1:28,0:13) - LOGICAL coffkind(1:28,0:13) - - COMPLEX, PARAMETER :: one = ( 1.0, 0.0) - COMPLEX, PARAMETER :: mone = (-1.0, 0.0) - COMPLEX, PARAMETER :: zro = ( 0.0, 0.0) - COMPLEX, PARAMETER :: ci = ( 0.0, 1.0) - COMPLEX, PARAMETER :: mci = ( 0.0,-1.0) - - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - ! ANGULAR MOMENTUM: 0 - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - - !*************************************************************************** - ! linear combinations of locs that are eigenfunctions - ! of z-reflection operation for l=0 (two atoms NOT in xy-plane) - !*************************************************************************** - ! even parity => leave as is - ! z-refl: even - DATA coffkind(1,0) /.true./ - DATA coeffi(1:2,1,0) /one,one/ - ! odd parity => multiply with i - ! z-refl: odd - DATA coffkind(2,0) /.false./ - DATA coeffi(1:2,2,0) /ci,mci/ - - - !*************************************************************************** - ! linear combinations of locs that are eigenfunctions - ! of z-reflection operation for l=0 (four atoms) - !*************************************************************************** - ! even parity => leave as is - ! z-refl: even - DATA coffkind(1,6) /.true./ - DATA coeffi(1:4,1,6) /one,one,one,one/ - - ! even parity => leave as is - ! z-refl: odd - DATA coffkind(2,6) /.false./ - DATA coeffi(1:4,2,6) /one,one,mone,mone/ - - ! odd parity => multiply with i - ! z-refl: even - DATA coffkind(3,6) /.true./ - DATA coeffi(1:4,3,6) /ci,mci,ci,mci/ - - ! odd parity => multiply with i - ! z-refl: odd - DATA coffkind(4,6) /.false./ - DATA coeffi(1:4,4,6) /ci,mci,mci,ci/ - - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - ! ANGULAR MOMENTUM: 1 - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - - - !******************************************************** - ! linear combinations of locs that are eigenfunctions - ! of z-reflection operation for l=1 (single atom) - !******************************************************** - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//odd parity => multiply with i - !z-refl: even - DATA coffkind(1,1) /.true./ - DATA coeffi(1:3,1,1) /one,zro,one/ - - !cos(theta) <==> Y10// odd parity => multiply with i - !z-refl: odd - DATA coffkind(2,1) /.false./ - DATA coeffi(1:3,2,1) /zro,ci,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// odd parity => multiply with i - !z-refl: even - DATA coffkind(3,1) /.true./ - DATA coeffi(1:3,3,1) /ci,zro,mci/ - - !******************************************************** - ! linear combinations of locs for l=1 (two atoms; both atoms in xy-plane) - !******************************************************** - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//even parity => leave as is - !z-refl: even - DATA coffkind(1,3) /.true./ - DATA coeffi(1:6,1,3) /ci,zro,ci,mci,zro,mci/ - - !cos(theta) <==> Y10//odd parity => multiply with i - !z-refl: odd - DATA coffkind(2,3) /.false./ - DATA coeffi(1:6,2,3) /zro,ci,zro,zro,ci,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// odd parity => multiply with i - !z-refl: even - DATA coffkind(3,3) /.true./ - DATA coeffi(1:6,3,3) /ci,zro,mci,ci,zro,mci/ - - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//odd parity => multiply with i - !z-refl: even - DATA coffkind(4,3) /.true./ - DATA coeffi(1:6,4,3) /one,zro,one,one,zro,one/ - - !cos(theta) <==> Y10//even parity => leave as is - !z-refl: odd - DATA coffkind(5,3) /.false./ - DATA coeffi(1:6,5,3) /zro,one,zro,zro,mone,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// even parity => leave as is - !z-refl: even - DATA coffkind(6,3) /.true./ - DATA coeffi(1:6,6,3) /one,zro,mone,mone,zro,one/ - - !******************************************************** - ! linear combinations of locs for l=1 (two atoms; atoms NOT in xy-plane) - !******************************************************** - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//even parity => leave as is - !z-refl: odd - DATA coffkind(1,4) /.false./ - DATA coeffi(1:6,1,4) /ci,zro,ci,mci,zro,mci/ - - !cos(theta) <==> Y10//odd parity => multiply with i - !z-refl: odd - DATA coffkind(2,4) /.false./ - DATA coeffi(1:6,2,4) /zro,ci,zro,zro,ci,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// odd parity => multiply with i - !z-refl: even - DATA coffkind(3,4) /.true./ - DATA coeffi(1:6,3,4) /ci,zro,mci,ci,zro,mci/ - - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//odd parity => multiply with i - !z-refl: even - DATA coffkind(4,4) /.true./ - DATA coeffi(1:6,4,4) /one,zro,one,one,zro,one/ - - !cos(theta) <==> Y10//even parity => leave as is - !z-refl: even - DATA coffkind(5,4) /.true./ - DATA coeffi(1:6,5,4) /zro,one,zro,zro,mone,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// even parity => leave as is - !z-refl: odd - DATA coffkind(6,4) /.false./ - DATA coeffi(1:6,6,4) /one,zro,mone,mone,zro,one/ - - - !******************************************************** - ! linear combinations of locs for l=1 (four atoms) - !******************************************************** - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//even parity => leave as is - !z-refl: even - DATA coffkind(1,7) /.false./ - DATA coeffi(1:12,1,7) /ci,zro,ci,mci,zro,mci, ci,zro,ci,mci,zro,mci/ - - !cos(theta) <==> Y10//odd parity => multiply with i - !z-refl: odd - DATA coffkind(2,7) /.false./ - DATA coeffi(1:12,2,7) /zro,ci,zro,zro,ci,zro, zro,ci,zro,zro,ci,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// odd parity => multiply with i - !z-refl: even - DATA coffkind(3,7) /.true./ - DATA coeffi(1:12,3,7) /ci,zro,mci,ci,zro,mci, ci,zro,mci,ci,zro,mci/ - - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//odd parity => multiply with i - !z-refl: odd - DATA coffkind(4,7) /.false./ - DATA coeffi(1:12,4,7) / one,zro, one, one,zro, one, mone,zro,mone,mone,zro,mone/ - - !cos(theta) <==> Y10//even parity => leave as is - !z-refl: odd - DATA coffkind(5,7) /.false./ - DATA coeffi(1:12,5,7) /zro, one,zro,zro,mone,zro, zro, one,zro,zro,mone,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// even parity => leave as is - !z-refl: even - DATA coffkind(6,7) /.true./ - DATA coeffi(1:12,6,7) /one,zro,mone,mone,zro,one, one,zro,mone,mone,zro,one/ - - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//even parity => leave as is - !z-refl: odd - DATA coffkind(7,7) /.false./ - DATA coeffi(1:12,7,7) / ci,zro, ci,mci,zro,mci, mci,zro,mci, ci,zro, ci/ - - !cos(theta) <==> Y10//odd parity => multiply with i - !z-refl: even - DATA coffkind(8,7) /.true./ - DATA coeffi(1:12,8,7) /zro, ci,zro,zro, ci,zro, zro,mci,zro,zro,mci,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// odd parity => multiply with i - !z-refl: odd - DATA coffkind(9,7) /.false./ - DATA coeffi(1:12,9,7) / ci,zro,mci, ci,zro,mci, mci,zro, ci,mci,zro, ci/ - - !sin(phi)sin(theta) <==> i(Y11 + Y1-1)//odd parity => multiply with i - !z-refl: true - DATA coffkind(10,7) /.true./ - DATA coeffi(1:12,10,7) /one,zro,one,one,zro,one, one,zro,one,one,zro,one/ - - !cos(theta) <==> Y10//even parity => leave as is - !z-refl: even - DATA coffkind(11,7) /.true./ - DATA coeffi(1:12,11,7) /zro, one,zro,zro,mone,zro, zro,mone,zro,zro, one,zro/ - - !cos(phi)sin(theta) <==> (Y11 - Y1-1)// even parity => leave as is - !z-refl: odd - DATA coffkind(12,7) /.false./ - DATA coeffi(1:12,12,7) / one,zro,mone,mone,zro, one, mone,zro, one, one,zro,mone/ - - - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - ! ANGULAR MOMENTUM: 2 - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - - - !******************************************************** - ! linear combinations of locs that are eigenfunctions - ! of z-reflection operation for l=2 (single atom) - !******************************************************** - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)// even parity => take as is - !z-refl: even - DATA coffkind(1,2) /.true./ - DATA coeffi(1:5,1,2) /one,zro,zro,zro,one/ - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(2,2) /.false./ - DATA coeffi(1:5,2,2) / zro ,one ,zro ,mone ,zro/ - - !(3*cos^2(theta)-1) <==> Y20// even parity => take as is - !z-refl: even - DATA coffkind(3,2) /.true./ - DATA coeffi(1:5,3,2) / zro ,zro ,one ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(4,2) /.false./ - DATA coeffi(1:5,4,2) / zro , ci ,zro , ci ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// even parity => take as is - !z-refl: even - DATA coffkind(5,2) /.true./ - DATA coeffi(1:5,5,2) / ci ,zro ,zro ,zro ,mci / - - !****************************************************************************** - ! linear combinations of locs that are eigenfunctions - ! of z-reflection operation for l=2; two atoms NOT in xy-plane - !****************************************************************************** - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)// even parity => take as is - !z-refl: even - DATA coffkind(1,8) /.true./ - DATA coeffi(1:10,1,8) / one ,zro ,zro ,zro ,one ,one ,zro ,zro ,zro ,one / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(2,8) /.false./ - DATA coeffi(1:10,2,8) / zro ,one ,zro ,mone ,zro ,zro ,one ,zro ,mone ,zro / - - !(3*cos^2(theta)-1) <==> Y20// even parity => take as is - !z-refl: even - DATA coffkind(3,8) /.true./ - DATA coeffi(1:10,3,8) / zro ,zro ,one ,zro ,zro ,zro ,zro ,one ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(4,8) /.false./ - DATA coeffi(1:10,4,8) / zro , ci ,zro , ci ,zro ,zro , ci ,zro , ci ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// even parity => take as is - !z-refl: even - DATA coffkind(5,8) /.true./ - DATA coeffi(1:10,5,8) / ci ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,mci / - - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)//odd parity => factor i - !z-refl: odd - DATA coffkind(6,8) /.false./ - DATA coeffi(1:10,6,8) / ci ,zro ,zro ,zro , ci ,mci ,zro ,zro ,zro ,mci / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// odd parity => factor i - !z-refl: even - DATA coffkind(7,8) /.true./ - DATA coeffi(1:10,7,8) / zro , ci ,zro ,mci ,zro ,zro ,mci ,zro , ci ,zro / - - !(3*cos^2(theta)-1) <==> Y20// odd parity => factor i - !z-refl: odd - DATA coffkind(8,8) /.false./ - DATA coeffi(1:10,8,8) / zro ,zro , ci ,zro ,zro ,zro ,zro ,mci ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// odd parity => factor i - !z-refl: even - DATA coffkind(9,8) /.true./ - DATA coeffi(1:10,9,8) / zro ,one ,zro ,one ,zro ,zro ,mone ,zro ,mone ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// odd parity => factor i - !z-refl: odd - DATA coffkind(10,8) /.false./ - DATA coeffi(1:10,10,8) / one ,zro ,zro ,zro ,mone ,mone ,zro ,zro ,zro ,one / - - !****************************************************************************** - ! linear combinations of locs that are eigenfunctions - ! of z-reflection operation for l=2; two atoms IN xy-plane - !****************************************************************************** - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)// odd parity => factor i - !z-refl: even - DATA coffkind(1,9) /.true./ - DATA coeffi(1:10,1,9) / ci ,zro ,zro ,zro , ci ,mci ,zro ,zro ,zro ,mci / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(2,9) /.false./ - DATA coeffi(1:10,2,9) / zro ,one ,zro ,mone ,zro ,zro ,one ,zro ,mone ,zro / - - !(3*cos^2(theta)-1) <==> Y20// even parity => take as is - !z-refl: even - DATA coffkind(3,9) /.true./ - DATA coeffi(1:10,3,9) / zro ,zro ,one ,zro ,zro ,zro ,zro ,one ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(4,9) /.false./ - DATA coeffi(1:10,4,9) / zro , ci ,zro , ci ,zro ,zro , ci ,zro , ci ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// even parity => take as is - !z-refl: even - DATA coffkind(5,9) /.true./ - DATA coeffi(1:10,5,9) / ci ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,mci / - - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)//even parity => leave as is - !z-refl: even - DATA coffkind(6,9) /.true./ - DATA coeffi(1:10,6,9) / one ,zro ,zro ,zro ,one ,one ,zro ,zro ,zro ,one / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// odd parity => factor i - !z-refl: odd - DATA coffkind(7,9) /.false./ - DATA coeffi(1:10,7,9) / zro , ci ,zro ,mci ,zro ,zro ,mci ,zro , ci ,zro / - - !(3*cos^2(theta)-1) <==> Y20// odd parity => factor i - !z-refl: even - DATA coffkind(8,9) /.true./ - DATA coeffi(1:10,8,9) / zro ,zro , ci ,zro ,zro ,zro ,zro ,mci ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// odd parity => factor i - !z-refl: odd - DATA coffkind(9,9) /.false./ - DATA coeffi(1:10,9,9) / zro ,one ,zro ,one ,zro ,zro ,mone ,zro ,mone ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// odd parity => factor i - !z-refl: even - DATA coffkind(10,9) /.true./ - DATA coeffi(1:10,10,9) / one ,zro ,zro ,zro ,mone ,mone ,zro ,zro ,zro ,one / - - - !****************************************************************************** - ! linear combinations of locs that are eigenfunctions - ! of z-reflection operation for l=2; four entangled atoms - !****************************************************************************** - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)// even parity => take as is - !z-refl: even - DATA coffkind(1,10) /.true./ - DATA coeffi(1:20,1,10) / one ,zro ,zro ,zro ,one ,one ,zro ,zro ,zro ,one& - ,one ,zro ,zro ,zro ,one ,one ,zro ,zro ,zro ,one / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(2,10) /.false./ - DATA coeffi(1:20,2,10) / zro ,one ,zro ,mone ,zro ,zro ,one ,zro ,mone ,zro& - ,zro ,one ,zro ,mone ,zro ,zro ,one ,zro ,mone ,zro / - - !(3*cos^2(theta)-1) <==> Y20// even parity => take as is - !z-refl: even - DATA coffkind(3,10) /.true./ - DATA coeffi(1:20,3,10) / zro ,zro ,one ,zro ,zro ,zro ,zro ,one ,zro ,zro& - ,zro ,zro ,one ,zro ,zro ,zro ,zro ,one ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// even parity => take as is - !z-refl: odd - DATA coffkind(4,10) /.false./ - DATA coeffi(1:20,4,10) / zro , ci ,zro , ci ,zro ,zro , ci ,zro , ci ,zro& - ,zro , ci ,zro , ci ,zro ,zro , ci ,zro , ci ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// even parity => take as is - !z-refl: even - DATA coffkind(5,10) /.true./ - DATA coeffi(1:20,5,10) / ci ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,mci& - , ci ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,mci / - - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)//odd parity => factor i - !z-refl: odd - DATA coffkind(6,10) /.false./ - DATA coeffi(1:20,6,10) / ci ,zro ,zro ,zro , ci ,mci ,zro ,zro ,zro ,mci& - ,mci ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro , ci / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// odd parity => factor i - !z-refl: odd - DATA coffkind(7,10) /.false./ - DATA coeffi(1:20,7,10) / zro , ci ,zro ,mci ,zro ,zro ,mci ,zro , ci ,zro& - ,zro , ci ,zro ,mci ,zro ,zro ,mci ,zro , ci ,zro / - - !(3*cos^2(theta)-1) <==> Y20// odd parity => factor i - !z-refl: odd - DATA coffkind(8,10) /.false./ - DATA coeffi(1:20,8,10) / zro ,zro , ci ,zro ,zro ,zro ,zro ,mci ,zro ,zro& - ,zro ,zro ,mci ,zro ,zro ,zro ,zro , ci ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// odd parity => factor i - !z-refl: even - DATA coffkind(9,10) /.true./ - DATA coeffi(1:20,9,10) / zro ,one ,zro ,one ,zro ,zro ,mone ,zro ,mone ,zro& - ,zro ,mone ,zro ,mone ,zro ,zro ,one ,zro ,one ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// odd parity => factor i - !z-refl: odd - DATA coffkind(10,10) /.false./ - DATA coeffi(1:20,10,10) / one ,zro ,zro ,zro ,mone ,mone ,zro ,zro ,zro ,one& - ,mone ,zro ,zro ,zro ,one ,one ,zro ,zro ,zro ,mone / - - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)// even parity => take as is - !z-refl: odd - DATA coffkind(11,10) /.false./ - DATA coeffi(1:20,11,10) / one ,zro ,zro ,zro ,one ,one ,zro ,zro ,zro ,one& - ,mone ,zro ,zro ,zro ,mone ,mone ,zro ,zro ,zro ,mone / - - !sin^2(theta)cos(2*phi) <==> (Y22 + Y2-2)//odd parity => factor i - !z-refl: odd - DATA coffkind(12,10) /.false./ - DATA coeffi(1:20,12,10) / ci ,zro ,zro ,zro , ci ,mci ,zro ,zro ,zro ,mci& - , ci ,zro ,zro ,zro , ci ,mci ,zro ,zro ,zro ,mci / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// even parity => take as is - !z-refl: even - DATA coffkind(13,10) /.true./ - DATA coeffi(1:20,13,10) / zro ,one ,zro ,mone ,zro ,zro ,one ,zro ,mone ,zro& - ,zro ,mone ,zro ,one ,zro ,zro ,mone ,zro ,one ,zro / - - !sin(theta)cos(theta)cos(phi) <==> (Y21 - Y2-1)// odd parity => factor i - !z-refl: even - DATA coffkind(14,10) /.true./ - DATA coeffi(1:20,14,10) / zro , ci ,zro ,mci ,zro ,zro ,mci ,zro , ci ,zro& - ,zro ,mci ,zro , ci ,zro ,zro , ci ,zro ,mci ,zro / - - !(3*cos^2(theta)-1) <==> Y20// even parity => take as is - !z-refl: odd - DATA coffkind(15,10) /.false./ - DATA coeffi(1:20,15,10) / zro ,zro ,one ,zro ,zro ,zro ,zro ,one ,zro ,zro& - ,zro ,zro ,mone ,zro ,zro ,zro ,zro ,mone ,zro ,zro / - - !(3*cos^2(theta)-1) <==> Y20// odd parity => factor i - !z-refl: even - DATA coffkind(16,10) /.true./ - DATA coeffi(1:20,16,10) / zro ,zro , ci ,zro ,zro ,zro ,zro ,mci ,zro ,zro& - ,zro ,zro , ci ,zro ,zro ,zro ,zro ,mci ,zro ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// even parity => take as is - !z-refl: even - DATA coffkind(17,10) /.true./ - DATA coeffi(1:20,17,10) / zro , ci ,zro , ci ,zro ,zro , ci ,zro , ci ,zro& - ,zro ,mci ,zro ,mci ,zro ,zro ,mci ,zro ,mci ,zro / - - !sin(theta)cos(theta)sin(phi) <==> i(Y21 + Y2-1)// odd parity => factor i - !z-refl: odd - DATA coffkind(18,10) /.false./ - DATA coeffi(1:20,18,10) / zro ,one ,zro ,one ,zro ,zro ,mone ,zro ,mone ,zro& - ,zro ,one ,zro ,one ,zro ,zro ,mone ,zro ,mone ,zro / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// even parity => take as is - !z-refl: odd - DATA coffkind(19,10) /.false./ - DATA coeffi(1:20,19,10) / ci ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,mci& - ,mci ,zro ,zro ,zro , ci ,mci ,zro ,zro ,zro , ci / - - !sin^2(theta)sin(2*phi) <==> i(Y22 - Y2-2)// odd parity => factor i - !z-refl: even - DATA coffkind(20,10) /.true./ - DATA coeffi(1:20,20,10) / one ,zro ,zro ,zro ,mone ,mone ,zro ,zro ,zro ,one& - ,one ,zro ,zro ,zro ,mone ,mone ,zro ,zro ,zro ,one / - - - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - ! ANGULAR MOMENTUM: 3 - !*************************************************************************** - !*************************************************************************** - !*************************************************************************** - - - !************************************************************* - ! linear combinations of locs that are eigenfunctions of - ! z-reflection operation for l=3 (single locs) - !************************************************************* - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(1,5) /.true./ - DATA coeffi(1:7,1,5) / one ,zro ,zro ,zro ,zro ,zro ,one / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(2,5) /.false./ - DATA coeffi(1:7,2,5) / zro , ci ,zro ,zro ,zro , ci ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(3,5) /.true./ - DATA coeffi(1:7,3,5) / zro ,zro ,one ,zro ,one ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//odd parity=>multiply by i - !z-refl: odd - DATA coffkind(4,5) /.false./ - DATA coeffi(1:7,4,5) / zro ,zro ,zro , ci ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(5,5) /.true./ - DATA coeffi(1:7,5,5) / zro ,zro , ci ,zro ,mci ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(6,5) /.false./ - DATA coeffi(1:7,6,5) / zro ,one ,zro ,zro ,zro ,mone ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(7,5) /.true./ - DATA coeffi(1:7,7,5) / ci ,zro ,zro ,zro ,zro ,zro ,mci / - - - - !*************************************************************************** - ! linear combinations of locs that are eigenfunctions of - ! z-reflection operation for l=3; two locs NOT in xy-plane (double locs) - !*************************************************************************** - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(1,11) /.true./ - DATA coeffi(1:14,1,11) / one ,zro ,zro ,zro ,zro ,zro ,one , one ,zro ,zro ,zro ,zro ,zro ,one / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(2,11) /.false./ - DATA coeffi(1:14,2,11) / zro , ci ,zro ,zro ,zro , ci ,zro , zro , ci ,zro ,zro ,zro , ci ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(3,11) /.true./ - DATA coeffi(1:14,3,11) / zro ,zro ,one ,zro ,one ,zro ,zro , zro ,zro ,one ,zro ,one ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//odd parity=>multiply by i - !z-refl: odd - DATA coffkind(4,11) /.false./ - DATA coeffi(1:14,4,11) / zro ,zro ,zro , ci ,zro ,zro ,zro , zro ,zro ,zro , ci ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(5,11) /.true./ - DATA coeffi(1:14,5,11) / zro ,zro , ci ,zro ,mci ,zro ,zro , zro ,zro , ci ,zro ,mci ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(6,11) /.false./ - DATA coeffi(1:14,6,11) / zro ,one ,zro ,zro ,zro ,mone ,zro , zro ,one ,zro ,zro ,zro ,mone ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(7,11) /.true./ - DATA coeffi(1:14,7,11) / ci ,zro ,zro ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,zro ,zro ,mci / - - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)//even parity => leave as is - !z-refl: odd - DATA coffkind(8,11) /.false./ - DATA coeffi(1:14,8,11) / ci ,zro ,zro ,zro ,zro ,zro , ci , mci ,zro ,zro ,zro ,zro ,zro ,mci / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)//even parity => leave as is - !z-refl: even - DATA coffkind(9,11) /.true./ - DATA coeffi(1:14,9,11) / zro ,one ,zro ,zro ,zro ,one ,zro , zro ,mone ,zro ,zro ,zro ,mone ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//even parity => leave as is - !z-refl: odd - DATA coffkind(10,11) /.false./ - DATA coeffi(1:14,10,11) / zro ,zro , ci ,zro , ci ,zro ,zro , zro ,zro ,mci ,zro ,mci ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//even parity => leave as is - !z-refl: even - DATA coffkind(11,11) /.true./ - DATA coeffi(1:14,11,11) / zro ,zro ,zro ,one ,zro ,zro ,zro , zro ,zro ,zro ,mone ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//even parity => leave as is - !z-refl: odd - DATA coffkind(12,11) /.false./ - DATA coeffi(1:14,12,11) / zro ,zro ,one ,zro ,mone ,zro ,zro , zro ,zro ,mone ,zro ,one ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)//even parity => leave as is - !z-refl: true - DATA coffkind(13,11) /.true./ - DATA coeffi(1:14,13,11) / zro , ci ,zro ,zro ,zro ,mci ,zro , zro ,mci ,zro ,zro ,zro , ci ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)//even parity => leave as is - !z-refl: odd - DATA coffkind(14,11) /.false./ - DATA coeffi(1:14,14,11) / one ,zro ,zro ,zro ,zro ,zro ,mone , mone ,zro ,zro ,zro ,zro ,zro ,one / - - !*************************************************************************** - ! linear combinations of locs that are eigenfunctions of - ! z-reflection operation for l=3; two locs IN xy-plane (double locs) - !*************************************************************************** - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(1,12) /.true./ - DATA coeffi(1:14,1,12) / one ,zro ,zro ,zro ,zro ,zro ,one , one ,zro ,zro ,zro ,zro ,zro ,one / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(2,12) /.false./ - DATA coeffi(1:14,2,12) / zro , ci ,zro ,zro ,zro , ci ,zro , zro , ci ,zro ,zro ,zro , ci ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(3,12) /.true./ - DATA coeffi(1:14,3,12) / zro ,zro ,one ,zro ,one ,zro ,zro , zro ,zro ,one ,zro ,one ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//odd parity=>multiply by i - !z-refl: odd - DATA coffkind(4,12) /.false./ - DATA coeffi(1:14,4,12) / zro ,zro ,zro , ci ,zro ,zro ,zro , zro ,zro ,zro , ci ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(5,12) /.true./ - DATA coeffi(1:14,5,12) / zro ,zro , ci ,zro ,mci ,zro ,zro , zro ,zro , ci ,zro ,mci ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(6,12) /.false./ - DATA coeffi(1:14,6,12) / zro ,one ,zro ,zro ,zro ,mone ,zro , zro ,one ,zro ,zro ,zro ,mone ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(7,12) /.true./ - DATA coeffi(1:14,7,12) / ci ,zro ,zro ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,zro ,zro ,mci / - - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)//even parity => leave as is - !z-refl: even - DATA coffkind(8,12) /.true./ - DATA coeffi(1:14,8,12) / ci ,zro ,zro ,zro ,zro ,zro , ci , mci ,zro ,zro ,zro ,zro ,zro ,mci / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)//even parity => leave as is - !z-refl: odd - DATA coffkind(9,12) /.false./ - DATA coeffi(1:14,9,12) / zro ,one ,zro ,zro ,zro ,one ,zro , zro ,mone ,zro ,zro ,zro ,mone ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//even parity => leave as is - !z-refl: even - DATA coffkind(10,12) /.true./ - DATA coeffi(1:14,10,12) / zro ,zro , ci ,zro , ci ,zro ,zro , zro ,zro ,mci ,zro ,mci ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//even parity => leave as is - !z-refl: odd - DATA coffkind(11,12) /.false./ - DATA coeffi(1:14,11,12) / zro ,zro ,zro ,one ,zro ,zro ,zro , zro ,zro ,zro ,mone ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//even parity => leave as is - !z-refl: even - DATA coffkind(12,12) /.true./ - DATA coeffi(1:14,12,12) / zro ,zro ,one ,zro ,mone ,zro ,zro , zro ,zro ,mone ,zro ,one ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)//even parity => leave as is - !z-refl: odd - DATA coffkind(13,12) /.false./ - DATA coeffi(1:14,13,12) / zro , ci ,zro ,zro ,zro ,mci ,zro , zro ,mci ,zro ,zro ,zro , ci ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)//even parity => leave as is - !z-refl: even - DATA coffkind(14,12) /.true./ - DATA coeffi(1:14,14,12) / one ,zro ,zro ,zro ,zro ,zro ,mone , mone ,zro ,zro ,zro ,zro ,zro ,one / - - !*************************************************************************** - ! linear combinations of locs that are eigenfunctions of - ! z-reflection operation for l=3; four entangled atoms - !*************************************************************************** - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(1,13) /.true./ - DATA coeffi(1:28,1,13) / one ,zro ,zro ,zro ,zro ,zro ,one , one ,zro ,zro ,zro ,zro ,zro ,one& - , one ,zro ,zro ,zro ,zro ,zro ,one , one ,zro ,zro ,zro ,zro ,zro ,one / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(2,13) /.false./ - DATA coeffi(1:28,2,13) / zro , ci ,zro ,zro ,zro , ci ,zro , zro , ci ,zro ,zro ,zro , ci ,zro& - , zro , ci ,zro ,zro ,zro , ci ,zro , zro , ci ,zro ,zro ,zro , ci ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(3,13) /.true./ - DATA coeffi(1:28,3,13) / zro ,zro ,one ,zro ,one ,zro ,zro , zro ,zro ,one ,zro ,one ,zro ,zro& - , zro ,zro ,one ,zro ,one ,zro ,zro , zro ,zro ,one ,zro ,one ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//odd parity=>multiply by i - !z-refl: odd - DATA coffkind(4,13) /.false./ - DATA coeffi(1:28,4,13) / zro ,zro ,zro , ci ,zro ,zro ,zro , zro ,zro ,zro , ci ,zro ,zro ,zro& - , zro ,zro ,zro , ci ,zro ,zro ,zro , zro ,zro ,zro , ci ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//odd parity=>multiply by i - !z-refl: even - DATA coffkind(5,13) /.true./ - DATA coeffi(1:28,5,13) / zro ,zro , ci ,zro ,mci ,zro ,zro , zro ,zro , ci ,zro ,mci ,zro ,zro& - , zro ,zro , ci ,zro ,mci ,zro ,zro , zro ,zro , ci ,zro ,mci ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(6,13) /.false./ - DATA coeffi(1:28,6,13) / zro ,one ,zro ,zro ,zro ,mone ,zro , zro ,one ,zro ,zro ,zro ,mone ,zro& - , zro ,one ,zro ,zro ,zro ,mone ,zro , zro ,one ,zro ,zro ,zro ,mone ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)// odd parity => multiply by i - !z-refl: even - DATA coffkind(7,13) /.true./ - DATA coeffi(1:28,7,13) / ci ,zro ,zro ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,zro ,zro ,mci& - , ci ,zro ,zro ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,zro ,zro ,mci / - - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)//even parity => leave as is - !z-refl: even - DATA coffkind(8,13) /.true./ - DATA coeffi(1:28,8,13) / ci ,zro ,zro ,zro ,zro ,zro , ci , mci ,zro ,zro ,zro ,zro ,zro ,mci& - , ci ,zro ,zro ,zro ,zro ,zro , ci , mci ,zro ,zro ,zro ,zro ,zro ,mci / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)//even parity => leave as is - !z-refl: odd - DATA coffkind(9,13) /.false./ - DATA coeffi(1:28,9,13) / zro ,one ,zro ,zro ,zro ,one ,zro , zro ,mone ,zro ,zro ,zro ,mone ,zro& - , zro ,one ,zro ,zro ,zro ,one ,zro , zro ,mone ,zro ,zro ,zro ,mone ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//even parity => leave as is - !z-refl: even - DATA coffkind(10,13) /.true./ - DATA coeffi(1:28,10,13) / zro ,zro , ci ,zro , ci ,zro ,zro , zro ,zro ,mci ,zro ,mci ,zro ,zro& - , zro ,zro , ci ,zro , ci ,zro ,zro , zro ,zro ,mci ,zro ,mci ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//even parity => leave as is - !z-refl: odd - DATA coffkind(11,13) /.false./ - DATA coeffi(1:28,11,13) / zro ,zro ,zro ,one ,zro ,zro ,zro , zro ,zro ,zro ,mone ,zro ,zro ,zro& - , zro ,zro ,zro ,one ,zro ,zro ,zro , zro ,zro ,zro ,mone ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//even parity => leave as is - !z-refl: even - DATA coffkind(12,13) /.true./ - DATA coeffi(1:28,12,13) / zro ,zro ,one ,zro ,mone ,zro ,zro , zro ,zro ,mone ,zro ,one ,zro ,zro& - , zro ,zro ,one ,zro ,mone ,zro ,zro , zro ,zro ,mone ,zro ,one ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)//even parity => leave as is - !z-refl: odd - DATA coffkind(13,13) /.false./ - DATA coeffi(1:28,13,13) / zro , ci ,zro ,zro ,zro ,mci ,zro , zro ,mci ,zro ,zro ,zro , ci ,zro& - , zro , ci ,zro ,zro ,zro ,mci ,zro , zro ,mci ,zro ,zro ,zro , ci ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)//even parity => leave as is - !z-refl: even - DATA coffkind(14,13) /.true./ - DATA coeffi(1:28,14,13) / one ,zro ,zro ,zro ,zro ,zro ,mone , mone ,zro ,zro ,zro ,zro ,zro ,one& - , one ,zro ,zro ,zro ,zro ,zro ,mone , mone ,zro ,zro ,zro ,zro ,zro ,one / - - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(15,13) /.false./ - DATA coeffi(1:28,15,13) / one ,zro ,zro ,zro ,zro ,zro ,one , one ,zro ,zro ,zro ,zro ,zro ,one& - , mone ,zro ,zro ,zro ,zro ,zro ,mone , mone ,zro ,zro ,zro ,zro ,zro ,mone / - - !sin^3(theta)sin(3*phi) <==> i(Y33+Y3-3)//even parity => leave as is - !z-refl: odd - DATA coffkind(16,13) /.false./ - DATA coeffi(1:28,16,13) / ci ,zro ,zro ,zro ,zro ,zro , ci , mci ,zro ,zro ,zro ,zro ,zro ,mci& - , mci ,zro ,zro ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,zro ,zro , ci / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)// odd parity => multiply by i - !z-refl: even - DATA coffkind(17,13) /.true./ - DATA coeffi(1:28,17,13) / zro , ci ,zro ,zro ,zro , ci ,zro , zro , ci ,zro ,zro ,zro , ci ,zro& - , zro ,mci ,zro ,zro ,zro ,mci ,zro , zro ,mci ,zro ,zro ,zro ,mci ,zro / - - !sin^2(theta)cos(theta)cos(2*phi)<==>(Y32+Y3-2)//even parity => leave as is - !z-refl: even - DATA coffkind(18,13) /.true./ - DATA coeffi(1:28,18,13) / zro ,one ,zro ,zro ,zro ,one ,zro , zro ,mone ,zro ,zro ,zro ,mone ,zro& - , zro ,mone ,zro ,zro ,zro ,mone ,zro , zro ,one ,zro ,zro ,zro ,one ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//odd parity=>multiply by i - !z-refl: odd - DATA coffkind(19,13) /.false./ - DATA coeffi(1:28,19,13) / zro ,zro ,one ,zro ,one ,zro ,zro , zro ,zro ,one ,zro ,one ,zro ,zro& - , zro ,zro ,mone ,zro ,mone ,zro ,zro , zro ,zro ,mone ,zro ,mone ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]sin(phi)<==>i(Y31+Y3-1)//even parity => leave as is - !z-refl: odd - DATA coffkind(20,13) /.false./ - DATA coeffi(1:28,20,13) / zro ,zro , ci ,zro , ci ,zro ,zro , zro ,zro ,mci ,zro ,mci ,zro ,zro& - , zro ,zro ,mci ,zro ,mci ,zro ,zro , zro ,zro , ci ,zro , ci ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//odd parity=>multiply by i - !z-refl: even - DATA coffkind(21,13) /.true./ - DATA coeffi(1:28,21,13) / zro ,zro ,zro , ci ,zro ,zro ,zro , zro ,zro ,zro , ci ,zro ,zro ,zro& - , zro ,zro ,zro ,mci ,zro ,zro ,zro , zro ,zro ,zro ,mci ,zro ,zro ,zro / - - ![5*cos^3(theta)-3*cos(theta)]<==>Y30//even parity => leave as is - !z-refl: even - DATA coffkind(22,13) /.true./ - DATA coeffi(1:28,22,13) / zro ,zro ,zro ,one ,zro ,zro ,zro , zro ,zro ,zro ,mone ,zro ,zro ,zro& - , zro ,zro ,zro ,mone ,zro ,zro ,zro , zro ,zro ,zro ,one ,zro ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//odd parity=>multiply by i - !z-refl: odd - DATA coffkind(23,13) /.false./ - DATA coeffi(1:28,23,13) / zro ,zro , ci ,zro ,mci ,zro ,zro , zro ,zro , ci ,zro ,mci ,zro ,zro& - , zro ,zro ,mci ,zro , ci ,zro ,zro , zro ,zro ,mci ,zro , ci ,zro ,zro / - - !sin(theta)[5*cos^2(theta)-1]cos(phi)<==>(Y31-Y3-1)//even parity => leave as is - !z-refl: odd - DATA coffkind(24,13) /.false./ - DATA coeffi(1:28,24,13) / zro ,zro ,one ,zro ,mone ,zro ,zro , zro ,zro ,mone ,zro ,one ,zro ,zro& - , zro ,zro ,mone ,zro ,one ,zro ,zro , zro ,zro ,one ,zro ,mone ,zro ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)// odd parity => multiply by i - !z-refl: even - DATA coffkind(25,13) /.true./ - DATA coeffi(1:28,25,13) / zro ,one ,zro ,zro ,zro ,mone ,zro , zro ,one ,zro ,zro ,zro ,mone ,zro& - , zro ,mone ,zro ,zro ,zro ,one ,zro , zro ,mone ,zro ,zro ,zro ,one ,zro / - - !sin^2(theta)cos(theta)sin(2*phi)<==>i(Y32-Y3-2)//even parity => leave as is - !z-refl: even - DATA coffkind(26,13) /.true./ - DATA coeffi(1:28,26,13) / zro , ci ,zro ,zro ,zro ,mci ,zro , zro ,mci ,zro ,zro ,zro , ci ,zro& - , zro ,mci ,zro ,zro ,zro , ci ,zro , zro , ci ,zro ,zro ,zro ,mci ,zro / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)// odd parity => multiply by i - !z-refl: odd - DATA coffkind(27,13) /.false./ - DATA coeffi(1:28,27,13) / ci ,zro ,zro ,zro ,zro ,zro ,mci , ci ,zro ,zro ,zro ,zro ,zro ,mci& - , mci ,zro ,zro ,zro ,zro ,zro , ci , mci ,zro ,zro ,zro ,zro ,zro , ci / - - !sin^3(theta)cos(3*phi) <==> (Y33 - Y3-3)//even parity => leave as is - !z-refl: odd - DATA coffkind(28,13) /.false./ - DATA coeffi(1:28,28,13) / one ,zro ,zro ,zro ,zro ,zro ,mone , mone ,zro ,zro ,zro ,zro ,zro ,one& - , mone ,zro ,zro ,zro ,zro ,zro ,one , one ,zro ,zro ,zro ,zro ,zro ,mone / - -END MODULE m_loccoeff - diff --git a/diagonalization/packed_to_full.f90 b/diagonalization/packed_to_full.f90 deleted file mode 100644 index 9e6ec09fb500d0ea79bede27efd4adef0704b607..0000000000000000000000000000000000000000 --- a/diagonalization/packed_to_full.f90 +++ /dev/null @@ -1,55 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Gruenberg Institut, Forschungszentrum Juelich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -MODULE m_packed_to_full -! Contains a service routine to expand a packed storage matrix to full storage -! A real and complex version is provided with a common interface - use m_juDFT - implicit none - private - interface packed_to_full - module procedure packed_to_full_r,packed_to_full_c - end interface packed_to_full - public packed_to_full -contains - subroutine packed_to_full_r(n,packed,full) - integer,intent(in) :: n - real,intent(in) :: packed(:) - real,allocatable,intent(out) :: full(:,:) - - integer:: i,err,i1,i2 - - ALLOCATE ( full(n,n), stat=err ) - if (err/=0) call judft_error("Allocation of full matrix failed",calledby="packed_to_full") - i=0 - DO i1=1,n - DO i2=1,i1 - i=i+1 - full(i2,i1)=packed(i) - full(i1,i2)=packed(i) - ENDDO - ENDDO - end subroutine packed_to_full_r - - subroutine packed_to_full_c(n,packed,full) - integer,intent(in) :: n - complex,intent(in) :: packed(:) - complex,allocatable,intent(out) :: full(:,:) - - integer:: i,err,i1,i2 - - ALLOCATE ( full(n,n), stat=err ) - if (err/=0) call judft_error("Allocation of full matrix failed",calledby="packed_to_full") - i=0 - DO i1=1,n - DO i2=1,i1 - i=i+1 - full(i2,i1)=packed(i) - full(i1,i2)=conjg(packed(i)) - ENDDO - ENDDO - end subroutine packed_to_full_c - end MODULE m_packed_to_full diff --git a/diagonalization/ssubredist1.F90 b/diagonalization/ssubredist1.F90 deleted file mode 100644 index 7a1432affa8ae30237245f0e97642a185ed75341..0000000000000000000000000000000000000000 --- a/diagonalization/ssubredist1.F90 +++ /dev/null @@ -1,671 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- -MODULE m_subredist1 - -CONTAINS -!DEC$ FREEFORM -SUBROUTINE subredist1(n,lda,SUB_COMM,nprow,npcol,iam,ierr,nb,achi_r,asca_r,achi_c,asca_c) - USE m_juDFT -#include"./cpp_double.h" - !#include"cpp_arch.h" - - IMPLICIT NONE - INCLUDE 'mpif.h' - - ! Will compute the redistribution of the half-colums processor iam owns - ! in the beginning when an n*n-matrix is distributed in the way - ! used in the beginning of chani: - ! Processor iam owns the upper half including the diagonal of columns - ! iam+(k*np)+1, k=0,...n/np[+1] - ! stored contigously in a one-dimensional array achi - ! to the matrix asca which is distributed to a - ! logical two-dimensional ScaLAPACK-grid with blocking size nb - ! For testing reasons nprow, npcol, n, and nb, the blocking size - ! for ScaLAPACK distribution - ! will be handed to the routine. Later on, nb will be computed from - ! nprow, npcol, and n. - ! The number of processors np=nprow*npcol will be computed - ! - ! Parameters, I/O channel numbers - REAL,OPTIONAL :: achi_r(:) ! Input, matrix in chani-distribution - REAL,OPTIONAL :: asca_r(:,:)!(lda,*) ! Output, matrix in ScaLAPACK distribution - COMPLEX,OPTIONAL :: achi_c(:) ! Input, matrix in chani-distribution - COMPLEX,OPTIONAL :: asca_c(:,:)!(lda,*) ! Output, matrix in ScaLAPACK distribution - ! Matrix might be real or complex - - INTEGER :: n,lda ! Global matrix size, local leading dimension of asca - INTEGER :: SUB_COMM,nprow,npcol,iam ! Communicator, - ! number of processor rows and columns in SUB_COMM, my rank - INTEGER :: nb ! blocking size, will be computed later - INTEGER :: ierr ! Error parameter to report problems - ! - ! Local Variables - INTEGER :: np ! number of procs, nprow*npcol=np - INTEGER :: me, myrow, mycol ! my processor coordinates - INTEGER :: j,i,k,l, kk ! counters - INTEGER :: mynumcols ! number of columns processor iam owns in the beginning - INTEGER :: jproccol, jprocrow ! processor column/row number counter variable - INTEGER :: isndproc, him ! processor in my column sending to me - INTEGER :: hisnumcols ! number of matrix columns proc isndproc owns - INTEGER :: iworldgroup ! mpi group index for group to SUB_COMM - INTEGER :: myrowssca,mycolssca ! number of rows/cols processor iam gets in Asca - INTEGER, ALLOCATABLE :: irowgroup(:), icolgroup(:) ! Subgroups rows, columns - INTEGER, ALLOCATABLE :: icommrow(:), icommcol(:) ! Communicators for Subgroups - INTEGER, ALLOCATABLE :: iranks(:) ! processor ranks to create subgroups - INTEGER, ALLOCATABLE :: glcolnum(:,:) ! Global column number of the columns I get - INTEGER :: jstart, jend, npfac ! see comments where used - INTEGER, ALLOCATABLE :: iwork(:,:), indx(:), isendcnts(:) ! see comments where used - INTEGER, ALLOCATABLE :: irecvcnts(:) ! see comments where used - INTEGER, ALLOCATABLE :: irecvcol(:) ! see comments where used - !Integer, allocatable :: ihavecolnum(:) ! Which columns do I have after first alltoall - INTEGER, ALLOCATABLE :: istartrecvbuf(:) ! Start of cols from proc i in recvbuf - INTEGER, ALLOCATABLE :: iinsendbuf(:) ! Now in sendbuf for proc i - INTEGER :: nsq ! Integer sqrt of np, to compute grid - INTEGER :: lensendbuf, lenrecvbuf ! Length of complete send- and recvbuf, first send - INTEGER :: lensendbuf2, lenrecvbuf2 ! Length of complete send- and recvbuf, second send - INTEGER :: jglob, jloc ! Global and local column index - INTEGER :: lglob, lloc ! Global and local row index - !Integer :: icntcol ! Number of columns I own after first alltoall - INTEGER :: indxrecvbuf ! Where I am in receivebuffer of first alltoall - INTEGER :: ierr1 ! For all the MPI calls - INTEGER :: subsize,mysubrank ! Test of comm_create - ! Arrays to redistribute - ! - REAL,ALLOCATABLE :: sendbuffer_r(:) - REAL,ALLOCATABLE :: recvbuffer_r(:) - COMPLEX,ALLOCATABLE :: sendbuffer_c(:) - COMPLEX,ALLOCATABLE :: recvbuffer_c(:) - - LOGICAL :: l_real - l_real=PRESENT(achi_r) - ! - ! Now do the computation - ! - np=nprow*npcol ! number of processors - mynumcols=n/np ! my number of columns in the chani-distribution - IF (iam < n-(n/np)*np) mynumcols=mynumcols+1 - myrow=iam/npcol ! my row number in the BLACS nprow*npcol grid - mycol=iam -(iam/npcol)*npcol ! my column number in the BLACS nprow*npcol grid - me = iam+1 ! makes life easyer - ! - ! Compute how many elements of Achi this processor owns (Never needed) - ! - ! j=0 - ! Do i=me,n,np ! Go through my columns of Achi - ! j=j+i ! Number of elements I have in Achi - ! Enddo - ! - ! Get Group handle for whole group - ! - CALL MPI_COMM_GROUP(SUB_COMM,iworldgroup, Ierr) - IF (ierr /= 0) THEN - WRITE(6,*) 'MPI_COMM_GROUP failed, ierr=',ierr - ! CALL CPP_flush(6) - CALL juDFT_error("MPI_COMM_GROUP failed") - ENDIF - ! - ! Allocate arrays to build subgroups for rows and columns - ! and create subgroups consisting of processor rows and columns - ! - ALLOCATE(iranks(1:MAX(npcol,nprow))) - ALLOCATE(irowgroup(0:nprow-1)) - ALLOCATE(icommrow(0:nprow-1)) - ALLOCATE(icolgroup(0:npcol-1)) - ALLOCATE(icommcol(0:npcol-1)) - DO k=0,nprow-1 ! Create groups of one processor row each - DO i=1,npcol - iranks(i)=k*npcol+i-1 - ENDDO - CALL MPI_GROUP_INCL(iworldgroup,npcol,iranks,irowgroup(k),& - ierr) - IF (ierr /= 0) THEN - WRITE(6,*) 'MPI_GROUP_INCL(',k,') row failed, ierr=',ierr - ! CALL CPP_flush(6) - CALL juDFT_error("MPI_COMM_INCL failed") - ENDIF - CALL MPI_COMM_CREATE(SUB_COMM,irowgroup(k),& - icommrow(k),ierr) - IF (ierr /= 0) THEN - WRITE(6,*) 'MPI_COMM_CREATE(',k,') row failed, ierr=',ierr - ! CALL CPP_flush(6) - CALL juDFT_error("MPI_COMM_CREATE failed") - ENDIF - ENDDO - DO k=0,npcol-1 ! Create groups of one processor column each - DO i=1,nprow - iranks(i)=k+(i-1)*npcol - ENDDO - CALL MPI_GROUP_INCL(iworldgroup,nprow,iranks,icolgroup(k),& - ierr) - IF (ierr /= 0) THEN - WRITE(6,*) 'MPI_GROUP_INCL(',k,') col failed, ierr=',ierr - !CALL CPP_flush(6) - CALL juDFT_error("MPI_COMM_INCL failed") - ENDIF - CALL MPI_COMM_CREATE(SUB_COMM,icolgroup(k),& - icommcol(k),ierr) - IF (ierr /= 0) THEN - WRITE(6,*) 'MPI_COMM_CREATE(',k,') col failed, ierr=',ierr - !CALL CPP_flush(6) - CALL juDFT_error("MPI_COMM_CREATE failed") - ENDIF - ENDDO - DEALLOCATE(iranks) - ! - ! - ! Allocate the arrays to store the indices - ! - ALLOCATE(indx(0:npcol-1)) - ALLOCATE(isendcnts(0:MAX(npcol-1,nprow-1))) - ALLOCATE(iinsendbuf(0:MAX(npcol-1,nprow-1))) - ! indx(jproccol) will indicate how many columns of processor column - ! jproccol are already found on this processor - indx=0 ! I have not found any column of any other processor yet - isendcnts=0 ! The length is 0 thus - ALLOCATE(iwork(1:2*mynumcols,0:npcol-1)) - iwork=0 - ! allocate(ihavecolnum( ((n/nb+1)/npcol+1)*nb )) - ! ihavecolnum=0 - ! iwork(2l-1,jproccol) will show the starting and iwork(2l,jproccol) - ! the end-indices for the columns of - ! processor column jproccol in the one-dimensional array of processor iam - ! Example: processor iam=me-1 has the upper parts of columns me+k*np, k=0,..,n/np[+1] - ! contigously in its local one-dimensional array. Thus elements - ! 1,..,me belong to global column jglob=me, elements me+1,...,2*me+np belong to - ! global column jglob=me+np ... - ! Global column jglob belongs to processor column - ! jproccol= modulo(jglob-1,nb*npcol)/nb - ! Thus iwork(1,modulo(me-1,nb*npcol)/nb)=1 - ! iwork(2,modulo(me-1,nb*npcol)/nb)=me - ! indx(modulo(me-1,nb*npcol)/nb)=indx(modulo(me-1,nb*npcol)/nb)+1 - ! iwork(2*indx(modulo(me+np-1,nb*npcol)/nb)+1,modulo(me+np-1,nb*npcol)/nb)=me+1 - ! iwork(2*indx(modulo(me+np-1,nb*npcol)/nb)+2,modulo(me+np-1,nb*npcol)/nb)=2*me+np - ! - jend=0 ! jstart will be jend+1 - npfac=0 ! my first elements belong to a column shorter or equal np - DO j=1,mynumcols ! Find out who my columns belong to - jstart=jend+1 ! Start point in achi - jend=jend+me+npfac*np ! End point in achi - jglob=me+npfac*np ! global column number - npfac=npfac+1 ! next column is np elements longer - jproccol=MODULO(jglob-1,npcol*nb)/nb ! Belongs to proc column - iwork(2*indx(jproccol)+1,jproccol)=jstart - iwork(2*indx(jproccol)+2,jproccol)=jend - isendcnts(jproccol)=isendcnts(jproccol)+jend-jstart+1 - ! Send isendcnts(jproccol) elements to (myrow,jproccol) - indx(jproccol)=indx(jproccol)+1 ! Send indx(jproccol) columns to him - ENDDO ! End j-loop over my columns - ! Count the number of columns and elements I have to send to which processor column - ! Here I only want to compute lensendbuf, thus I don't need to store - ! anything to sendbuffer because sendbuffer is not yet allocated!!! - lensendbuf=0 ! Length of sendbuffer - DO j=0,npcol-1 - IF (j /= mycol) lensendbuf=lensendbuf+isendcnts(j) - ENDDO - IF (l_real) THEN - ALLOCATE(sendbuffer_r(lensendbuf)) - ELSE - ALLOCATE(sendbuffer_c(lensendbuf)) - ENDIF - - ! - ! Now I have to put the elements to be sent to sendbuffer !! - ! - lensendbuf=0 - iinsendbuf(0)=0 ! Starting addresses in sendbuffer (See MPI_ALLTOALLV) - DO j=0,npcol-1 - IF (indx(j)==0) THEN - IF (isendcnts(j) /= 0) WRITE(6,*) iam, & - ' Fehler, isendcnts nicht 0',isendcnts(j) - IF (j < npcol-1) iinsendbuf(j+1)=iinsendbuf(j) - ELSE - IF (j /= mycol) THEN - IF (l_real) THEN - DO i=1,indx(j) - DO kk=iwork(2*i-1,j),iwork(2*i,j) - sendbuffer_r(lensendbuf+1+kk-iwork(2*i-1,j))=& - Achi_r(kk) - ENDDO - lensendbuf=lensendbuf+iwork(2*i,j)-iwork(2*i-1,j)+1 - ENDDO - ELSE - DO i=1,indx(j) - DO kk=iwork(2*i-1,j),iwork(2*i,j) - sendbuffer_c(lensendbuf+1+kk-iwork(2*i-1,j))=& - Achi_c(kk) - ENDDO - lensendbuf=lensendbuf+iwork(2*i,j)-iwork(2*i-1,j)+1 - ENDDO - ENDIF - ENDIF - IF (j == mycol) isendcnts(j)=0 ! must be 0 for mpi_alltoallv - IF (j < npcol-1) iinsendbuf(j+1)=iinsendbuf(j)+isendcnts(j) - ENDIF - ENDDO - ALLOCATE(irecvcnts(0:MAX(npcol-1,nprow-1))) ! Number of elements I will get - ALLOCATE(irecvcol(0:npcol-1)) ! Number of columns I will get - ALLOCATE(glcolnum(n/np+1,0:npcol-1)) ! Global column numbers - ALLOCATE(istartrecvbuf(0:MAX(npcol-1,nprow-1))) ! Start in receivebuffer - irecvcnts=0 - irecvcol=0 - glcolnum=0 - lenrecvbuf=0 - istartrecvbuf(0)=0 - DO k=0,npcol-1 ! Now look what processors in my row send to me - isndproc=myrow*npcol+k ! Number of processor (myrow,k) - IF (k /= mycol) THEN ! All except me are interesting - him=isndproc+1 ! make life easier - hisnumcols=n/np ! his number of columns in the chani-distribution - IF (isndproc < n-(n/np)*np) hisnumcols=hisnumcols+1 - jend=0 ! jstart will be jend+1 - npfac=0 ! his first elements belong to a column shorter or equal np - DO j=1,hisnumcols ! Look through his columns - jstart=jend+1 ! Starting in his part of Achi - jend=jend+him+npfac*np ! Ending in his part of Achi - jglob=him+npfac*np ! Global column number - npfac=npfac+1 ! next column is np elements longer - jproccol=MODULO(jglob-1,npcol*nb)/nb ! Will belong to proc column - IF (jproccol==mycol) THEN ! I will get this column - irecvcol(k)=irecvcol(k)+1 ! I get one more column from (myrow,k) - irecvcnts(k)=irecvcnts(k)+jend-jstart+1 - glcolnum(irecvcol(k),k)=jglob ! It is global column number - ENDIF ! I get this column - ENDDO ! j=1,hisnumcols - lenrecvbuf=lenrecvbuf+irecvcnts(k) - IF (k < npcol-1) istartrecvbuf(k+1) = & - istartrecvbuf(k)+irecvcnts(k) - ! what I got from k starts at istartrecvbuf(k)+1 in the recvbuf - ELSE ! it is me, start in recvbuf for the next is same as for this - IF (k < npcol-1) istartrecvbuf(k+1)= istartrecvbuf(k) - ENDIF ! If it is not me - ENDDO ! Loop over procs in my row - IF (l_real) THEN - ALLOCATE(recvbuffer_r(lenrecvbuf)) - CALL MPI_ALLTOALLV(sendbuffer_r,isendcnts,iinsendbuf,& - CPP_MPI_TYP_REAL,recvbuffer_r,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_REAL,icommrow(myrow),ierr) - ELSE - ALLOCATE(recvbuffer_c(lenrecvbuf)) - CALL MPI_ALLTOALLV(sendbuffer_c,isendcnts,iinsendbuf,& - CPP_MPI_TYP_COMPLEX,recvbuffer_c,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_COMPLEX,icommrow(myrow),ierr) - ENDIF - ! - ! Now look to Asca to put the elements of Achi or receivebuffer to - ! See later, where the number of elements of column jglob this - ! processor will get is computed. Jglob there plays the same role as - ! n here. - ! - myrowssca=(n-1)/(nb*nprow)*nb+ & - MIN(MAX(n-(n-1)/(nb*nprow)*nb*nprow-nb*myrow,0),nb) - ! My number of rows in ScaLAPACK distribution - IF (myrowssca > lda) & - WRITE(6,*) 'Redist1: Wrong dimension of Asca, lda=',lda, & - ' myrowssca=',myrowssca - ! My number of columns in ScaLAPACK distribution - mycolssca=(n-1)/(nb*npcol)*nb+ & - MIN(MAX(n-(n-1)/(nb*npcol)*nb*npcol-nb*mycol,0),nb) - ! - ! See computing jloc from jglob - ! - istartrecvbuf=istartrecvbuf+1 ! Makes live easier - ! - ! Now look what this processor should have after the first alltoall - ! along its processor row - ! - ! First those columns it already had, they are still on Achi - ! - isendcnts=0 ! Number of elements alltoall to send along column - ! icntcol=0 - lensendbuf2=0 ! New sendbuffer - DO k=1,indx(mycol) ! I had indx(mycol) columns for my processor column - ! I have to store it to Asca or sent in next Alltoallv - jglob = iwork(2*k,mycol)-iwork(2*k-1,mycol)+1 !global column index - jloc=jglob-(jglob-1)/(nb*npcol)*nb*(npcol-1)-nb*mycol - ! jloc is local column index for ScaLAPACK distribution - ! Achtung, hier stand urspruenglich - ! jloc= jglob - jglob/(nb*npcol)*nb*(npcol-1)-nb*mycol, aber das - ! fuehrte im Fall der letzten Spalte eines "sweeps" zu einem Fehler - ! Gedanke, der diese Formel ergibt: - ! Man treagt von der globalen Spaltennummer jeweils einen kompletten - ! "sweep" von nb*npcol Spalten ab, also hat man - ! jglob/(nb*npcol) mal jedem Prozessor dieser Zeile einen Spalten - ! block gegeben, also habe auch ich lokal jglob/(nb*npcol) ganze - ! Spaltenbloecke, es beleiben - ! jglob - jglob/(nb*npcol)*nb*npcol Spalten uebrig, die noch - ! zu verteilen sind, davon haben die mycol Prozessoren, die vor - ! mir in der Reihe sind jeweils einen ganzen Block, - ! bleiben also jglob - jglob/(nb*npcol)*nb*npcol - nb*mycol - ! Elemente uebrig. - ! Damit bin ich lokal bei folgender Spalte: - ! jglob/(nb*npcol) ganze "sweeps", also lokale Bloecke bis zu dieser - ! Spalte, das sind also jglob/(nb*npcol)*nb Spalten, dann noch die Restspalten, - ! Also jglob - jglob/(nb*npcol)*nb*npcol - nb*mycol Spalten. - ! Das geht gut, ausser bei der letzten Spalte eines "sweeps". - ! Hier ergibt sich jglob - jglob/(nb*npcol)*nb*npcol=0, also keine Spalten - ! uebrig. Diese letzte Spalte hat aber Prozessor (myrow, npcol-1), also - ! mycol=npcol-1, das ist meistens groesser als 0, also werden - ! nb*(npcol-1) Spalten vom Rest abgezogen, dieser Rest ist aber =0, also - ! entsteht eine negative Zahl, die zu den Spalten von den ganzen Bloecken - ! addiert wird. Eigentlich darf aber nichts mehr addiert oder subtrahiert - ! Werden, ich habe genau jglob/(nb*npcol)*nb Spalten. - ! Daher muss die Zahl der kompletten "sweeps" berechnet werden, die vor - ! dieser Spalte jglob schon fertig verteilt sind, dies ist dieselbe Anzahl - ! ausser in dem einen Fall, in dem es gerade aufgeht. - ! Wenn man dann jglob - (jglob-1)/(nb*npcol)*nb*npcol berechnet, kommt - ! genau nb*npcol, also ein ganzer "sweep" als Rest heraus. - ! Zieht man jetzt die Spalten ab, die die Prozessoren vor mir, - ! also vor (myrow, npcol-1) haben, so erhaelt man - ! nb*npcol-nb*mycol=nb*(npcol-mycol)=nb*(npcol-(npcol-1))=nb als lokalen Index - ! fuer den Rest, und ich bekomme als lokalen Spaltenindex - ! (Anzahl kompletter "sweeps" bis jglob-1)*nb+nb - ! =(jglob-1)/(nb*npcol)*nb+nb, und da in diesem speziellen Fall, und nur - ! in diesem Fall gilt (jglob-1)/(nb*npcol)=jglob/(nb*npcol)-1 - ! erhalte ich jloc=jglob/(nb*npcol)*nb, und das stimmt, denn es ist die - ! letzte Spalte im Block jglob/(nb*npcol) - lloc=1 ! local row index for ScaLAPACK distribution - ! icntcol= icntcol+1 - ! ihavecolnum(icntcol)=jglob - DO j=1,jglob,nb - ! Go through column jglob by row blocks - ! The row block treated now belongs to processor row - ! mod(j/nb,nprow) - IF (MOD(j/nb,nprow)==myrow) THEN ! Belonging to my processor row - ! = I get it, store it to Asca - IF (l_real) THEN - DO lglob=j,MIN(jglob,j+nb-1) ! This row block, no more than jglob - Asca_r(lloc,jloc)=Achi_r(iwork(2*k-1,mycol)+lglob-1) - lloc=lloc+1 - ENDDO - ELSE - DO lglob=j,MIN(jglob,j+nb-1) ! This row block, no more than jglob - Asca_c(lloc,jloc)=Achi_c(iwork(2*k-1,mycol)+lglob-1) - lloc=lloc+1 - ENDDO - END IF - ELSE ! send to processor row mod(j/nb,nprow) - isendcnts(MOD(j/nb,nprow))=isendcnts(MOD(j/nb,nprow))& - + MIN(nb,jglob-j+1) ! Compute only number of elements - ENDIF - ENDDO - ENDDO - ! - ! Now processing the received columns, processor by processor - ! - ! icntcol=indx(mycol) - indxrecvbuf=0 ! Go through receivebuffer element by element - DO l=0,npcol-1 - DO k=1,irecvcol(l) ! If l=mycol irecvcol=0, nothing will be done - jglob=glcolnum(k,l) ! Received column is global column number jglob - jloc=jglob-(jglob-1)/(nb*npcol)*nb*(npcol-1)-nb*mycol - ! It will be local column number jloc of Asca - ! icntcol = icntcol+1 ! Number of columns I have now - ! ihavecolnum(icntcol)=jglob ! Global indices of columns I have - lloc=1 ! Now go through column by row blocks - DO j=1,jglob,nb - ! Go through column jglob by row blocks - ! The row block treated now belongs to processor row - ! mod(j/nb,nprow) - IF (MOD(j/nb,nprow)==myrow) THEN ! my row block, store it to Asca - IF (l_real) THEN - DO lglob=j,MIN(jglob,j+nb-1) ! This row block, no more than jglob - indxrecvbuf=indxrecvbuf+1 - Asca_r(lloc,jloc)=recvbuffer_r(indxrecvbuf) - lloc=lloc+1 - ENDDO - ELSE - DO lglob=j,MIN(jglob,j+nb-1) ! This row block, no more than jglob - indxrecvbuf=indxrecvbuf+1 - Asca_c(lloc,jloc)=recvbuffer_c(indxrecvbuf) - lloc=lloc+1 - ENDDO - ENDIF - ELSE ! not my row block send to processor row mod(j/nb,nprow) - ! min(nb,jglob-j+1) more elements - isendcnts(MOD(j/nb,nprow))= & - isendcnts(MOD(j/nb,nprow))+ MIN(nb,jglob-j+1) - indxrecvbuf=indxrecvbuf+MIN(nb,jglob-j+1) - ENDIF ! who's row block - ENDDO ! End going through column jglob - ENDDO ! End going through columns received from processor myrow,k - ENDDO ! End going through all processors in my row - ! Deallocate some workspace no longer needed - DEALLOCATE(glcolnum) - DEALLOCATE(irecvcol) - ! deallocate(ihavecolnum) - ! - iinsendbuf(0)=1 - DO k=0,nprow-1 - IF (k /= myrow) THEN ! I will send isendcnts(k )elements to processor - ! k,mycol - lensendbuf2=lensendbuf2+isendcnts(k) - IF (k < nprow-1) iinsendbuf(k+1)=iinsendbuf(k)+isendcnts(k) - ELSE ! It is my row, don't have to send - IF (k < nprow-1) iinsendbuf(k+1)=iinsendbuf(k) - ENDIF - ENDDO - IF (lensendbuf2 > lensendbuf) THEN - IF (l_real) THEN - DEALLOCATE(sendbuffer_r) - ALLOCATE(sendbuffer_r(lensendbuf2)) - ELSE - DEALLOCATE(sendbuffer_c) - ALLOCATE(sendbuffer_c(lensendbuf2)) - ENDIF - ENDIF - ! - ! Now look through all columns I need and see whether I have them or get them - ! - irecvcnts=0 ! Count for next receive - lenrecvbuf2=0 ! New receivebuffer - jstart=0 ! No column from Achi stored to sendbuf - DO k=mycol*nb+1,n,npcol*nb ! k is the first element in a column block - DO jglob=k,MIN(k+nb-1,n) ! global column index - ! the processor originally having this column is - ! processor modulo(jglob-1,np) - IF (MODULO(jglob-1,np)==iam) THEN ! It was my column in the beginning - jstart=jstart+1 - DO j=1,jglob,nb ! Row block starting with j - ! Go through column jglob by row blocks - ! The row block treated now belongs to processor row - ! mod(j/nb,nprow) - IF (MOD(j/nb,nprow) /= myrow) THEN ! not my row block - ! Now find out where it is in Achi - IF (l_real) THEN - DO l=j,MIN(j+nb-1,jglob) - sendbuffer_r(iinsendbuf(MOD(j/nb,nprow)))=& - Achi_r(iwork(2*jstart-1,mycol)+l-1) - iinsendbuf(MOD(j/nb,nprow))= & - iinsendbuf(MOD(j/nb,nprow))+1 - ENDDO - ELSE - DO l=j,MIN(j+nb-1,jglob) - sendbuffer_c(iinsendbuf(MOD(j/nb,nprow)))=& - Achi_c(iwork(2*jstart-1,mycol)+l-1) - iinsendbuf(MOD(j/nb,nprow))= & - iinsendbuf(MOD(j/nb,nprow))+1 - ENDDO - ENDIF - ENDIF ! who's row block - ENDDO ! End going through column jglob - ELSE IF (MODULO(jglob-1,np)/npcol==myrow) THEN - jproccol=MODULO(MODULO(jglob-1,np),npcol) - ! modulo(jglob-1,np) is the number of the processor who had - ! this column in the beginning - ! This processor is in my processor row, if its number divided by - ! npcol is the same as my number divided by npcol, and this is - ! myrow - ! processor number divided by npcol = processor row number - ! processor number modulo npcol = processor column number - ! It is in my receive buffer - DO j=1,jglob,nb - ! Go through column jglob by row blocks - ! The row block treated now belongs to processor row - ! mod(j/nb,nprow) - IF (MOD(j/nb,nprow)==myrow) THEN ! my row block - ! stored the elements from recvbuf to Asca - istartrecvbuf(jproccol)=istartrecvbuf(jproccol)+& - MIN(nb,jglob-j+1) - ELSE ! not my row block, store to sendbuffer for row mod(j/nb,nprow) - ! Now find out where it is in recvbuffer - DO l=j,MIN(j+nb-1,jglob) - IF (l_real) THEN - sendbuffer_r(iinsendbuf(MOD(j/nb,nprow)))=recvbuffer_r(istartrecvbuf(jproccol)) - ELSE - sendbuffer_c(iinsendbuf(MOD(j/nb,nprow)))=recvbuffer_c(istartrecvbuf(jproccol)) - END IF - - iinsendbuf(MOD(j/nb,nprow))= & - iinsendbuf(MOD(j/nb,nprow))+1 - istartrecvbuf(jproccol)= & - istartrecvbuf(jproccol)+1 - ENDDO - ENDIF ! who's row block - ENDDO ! End going through column jglob - ELSE ! I get it in this alltoall - jprocrow=MODULO(jglob-1,np)/npcol ! processor row where it is - ! The number of elements I have to receive for this column is: - ! The number of complete "sweeps" of nb*nprow elements of this - ! column are (jglob-1)/(nb*nprow), that is, each processor in my - ! processor column gets at least (jglob-1)/(nb*nprow)*nb elements - ! of column jglob. - ! The rest is jglob-(jglob-1)/(nb*nprow)*nb*nprow elements. - ! Out of these elements I will get - ! at least 0, at most nb and, if there is just part of a block - ! for me, I'll get rest-myrow*nb elements more - ! Thus the elements I get in addition to the (jglob-1)/(nb*nprow)*nb - ! each processor gets is - ! min(max(0,jglob-(jglob-1)/(nb*nprow)*nb*nprow-myrow*nb),nb) - irecvcnts(jprocrow)=irecvcnts(jprocrow)+& - (jglob-1)/(nb*nprow)*nb + & - MIN(MAX(0,jglob-(jglob-1)/(nb*nprow)*nb*nprow-myrow*nb),nb) - ENDIF - ENDDO ! jglob=k,min(k+nb-1,n) - ENDDO ! k=mycol*nb+1,n,npcol*nb - ! Compute length of new receivebuffer - ! and start addresses in the receivebuffer - istartrecvbuf(0)=0 - iinsendbuf(0)=0 - DO jprocrow=0,nprow-1 - IF (jprocrow /= myrow) THEN ! I will get and send in this alltoall - lenrecvbuf2=lenrecvbuf2+irecvcnts(jprocrow) - IF (jprocrow < nprow-1) THEN - istartrecvbuf(jprocrow+1)=istartrecvbuf(jprocrow)+& - irecvcnts(jprocrow) - iinsendbuf(jprocrow+1)=iinsendbuf(jprocrow)+ & - isendcnts(jprocrow) - ENDIF - ELSE ! this is my row, i won't get anything but have to count - ! start in receivebuffer - IF (jprocrow < nprow-1) THEN - istartrecvbuf(jprocrow+1)= istartrecvbuf(jprocrow) - iinsendbuf(jprocrow+1)=iinsendbuf(jprocrow) - ENDIF - ENDIF - ENDDO - IF (lenrecvbuf2 > lenrecvbuf) THEN - IF (l_real) THEN - DEALLOCATE(recvbuffer_r) ! Old one no longer needed if it is too small - ALLOCATE(recvbuffer_r(lenrecvbuf2)) - ELSE - DEALLOCATE(recvbuffer_c) ! Old one no longer needed if it is too small - ALLOCATE(recvbuffer_c(lenrecvbuf2)) - ENDIF - ENDIF - IF (l_real) THEN - CALL MPI_ALLTOALLV(sendbuffer_r,isendcnts,iinsendbuf,& - CPP_MPI_TYP_REAL,recvbuffer_r,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_REAL,icommcol(mycol),ierr) - ELSE - CALL MPI_ALLTOALLV(sendbuffer_c,isendcnts,iinsendbuf,& - CPP_MPI_TYP_COMPLEX,recvbuffer_c,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_COMPLEX,icommcol(mycol),ierr) - ENDIF - IF (l_real) THEN - DEALLOCATE(sendbuffer_r) - ELSE - DEALLOCATE(sendbuffer_c) - END IF - - istartrecvbuf=istartrecvbuf+1 - ! - ! After second alltoall store elements received to Asca - ! Go again through colums I have to get - ! - DO k=mycol*nb+1,n,npcol*nb ! k is the first element in a column block - DO jglob=k,MIN(k+nb-1,n) ! global column index - ! the processor originally having this column is - ! processor modulo(jglob-1,np) - IF (MODULO(jglob-1,np)/npcol /= myrow) THEN ! I received it in second alltoallv - ! modulo(jglob-1,np) is the number of the processor who had - ! this column in the beginning - ! This processor is in my processor row, if its number divided by - ! npcol is the same as my number divided by npcol, and this is - ! myrow - ! In that case nothing has to be done, I already treated this - ! column before second alltoall - jprocrow=MODULO(jglob-1,np)/npcol - ! The number of elements I have to receive for this column is: - ! The number of complete "sweeps" of nb*nprow elements of this - ! column are (jglob-1)/(nb*nprow), that is, each processor in my - ! processor column gets at least (jglob-1)/(nb*nprow)*nb elements - ! of column jglob. - ! The rest is jglob-(jglob-1)/(nb*nprow)*nb*nprow elements. - ! Out of these elements I will get - ! at least 0, at most nb and, if there is just part of a block - ! for me, I'll get rest-myrow*nb elements more - ! Thus the elements I get in addition to the (jglob-1)/(nb*nprow)*nb - ! each processor gets is - ! min(max(0,jglob-(jglob-1)/(nb*nprow)*nb*nprow-myrow*nb),nb) - jloc=jglob-(jglob-1)/(nb*npcol)*nb*(npcol-1)-nb*mycol - IF (l_real) THEN - DO i=1,(jglob-1)/(nb*nprow)*nb + & - MIN(MAX(0,jglob-(jglob-1)/(nb*nprow)*nb*nprow-myrow*nb),nb) - Asca_r(i,jloc)=recvbuffer_r(istartrecvbuf(jprocrow)+i-1) - ENDDO - ELSE - DO i=1,(jglob-1)/(nb*nprow)*nb + & - MIN(MAX(0,jglob-(jglob-1)/(nb*nprow)*nb*nprow-myrow*nb),nb) - Asca_c(i,jloc)=recvbuffer_c(istartrecvbuf(jprocrow)+i-1) - ENDDO - END IF - istartrecvbuf(jprocrow) = istartrecvbuf(jprocrow)+& - (jglob-1)/(nb*nprow)*nb + & - MIN(MAX(0,jglob-(jglob-1)/(nb*nprow)*nb*nprow-myrow*nb),nb) - ENDIF ! End if I got it in second alltoallv - ENDDO ! jglob=k,min(k+nb-1,n) - ENDDO ! k=mycol*nb+1,n,npcol*nb - ! - ! Free the row and column communicators - ! - DO k=0,nprow-1 - IF (icommrow(k) /= MPI_COMM_NULL) THEN - CALL MPI_COMM_FREE(icommrow(k),ierr) - ENDIF - CALL MPI_GROUP_FREE(irowgroup(k),ierr) - ENDDO - DO k=0,npcol-1 - IF (icommcol(k) /= MPI_COMM_NULL) THEN - CALL MPI_COMM_FREE(icommcol(k),ierr) - ENDIF - CALL MPI_GROUP_FREE(icolgroup(k),ierr) - ENDDO - ! - ! Deallocate the arrays to store the indices - ! - DEALLOCATE(irecvcnts) - DEALLOCATE(istartrecvbuf) - DEALLOCATE(iinsendbuf) - DEALLOCATE(indx) - DEALLOCATE(isendcnts) - DEALLOCATE(iwork) - DEALLOCATE(irowgroup) - DEALLOCATE(icommrow) - DEALLOCATE(icolgroup) - DEALLOCATE(icommcol) - RETURN -END SUBROUTINE subredist1 -END diff --git a/diagonalization/ssubredist2.F90 b/diagonalization/ssubredist2.F90 deleted file mode 100644 index 1d7498b4f9240d14ce5635d66b0ea4138aa929c0..0000000000000000000000000000000000000000 --- a/diagonalization/ssubredist2.F90 +++ /dev/null @@ -1,521 +0,0 @@ -MODULE m_subredist2 -CONTAINS -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -!DEC$ FREEFORM -SUBROUTINE subredist2(n,m,lda,SUB_COMM,nprow,npcol,& - iam,ierr,nb,achi_r,asca_r,achi_c,asca_c) -#include"./cpp_double.h" - !#include"cpp_arch.h" - USE m_juDFT - IMPLICIT NONE - INCLUDE 'mpif.h' - ! Will compute the redistribution of the n*m-matrix - ! Asca distributed in a block-cyclic 2-dimensional manner with - ! block-size nb coming from ScaLAPACK - ! to a column cyclically distributed matrix - ! Achi used as output of chani. - ! nprow, npcol, n,m, and nb, the blocking size for ScaLAPACK distribution - ! will be handed to the routine. - ! The number of processors np=nprow*npcol will be computed - ! - ! Parameters, I/O channel numbers - REAL,OPTIONAL :: achi_r(*) ! Output, matrix in row cyclic-distribution - REAL,OPTIONAL :: asca_r(lda,*) ! Input, matrix in ScaLAPACK distribution - COMPLEX,OPTIONAL :: achi_c(*) ! Input, matrix in chani-distribution - COMPLEX,OPTIONAL :: asca_c(lda,*) ! Output, matrix in ScaLAPACK distribution - - INTEGER :: n,m,lda ! Global matrix sizes, local leading dimension of asca - INTEGER :: SUB_COMM,nprow,npcol,iam ! Communicator, - ! number of processor rows and columns in SUB_COMM, my rank - INTEGER :: nb ! blocking size, will be computed later - INTEGER :: ierr ! Error parameter to report problems - ! - ! Local Variables - INTEGER :: np ! number of procs, nprow*npcol=np - INTEGER :: myrow, mycol ! my processor coordinates - INTEGER :: j,i,k,l, kk ! counters - INTEGER :: mynumcols ! number of columns processor iam owns after redistribution - INTEGER :: jproc, jprocrow ! processor (row) number counter variable - INTEGER :: hiscolssca,hisrowssca ! number of matrix columns/rows sending proc - INTEGER :: iworldgroup ! mpi group index for group to SUB_COMM - INTEGER :: myrowssca,mycolssca ! number of rows/cols processor iam has in Asca - INTEGER, ALLOCATABLE :: irowgroup(:), icolgroup(:) ! Subgroups rows, columns - INTEGER, ALLOCATABLE :: icommrow(:), icommcol(:) ! Communicators for Subgroups - INTEGER, ALLOCATABLE :: iranks(:) ! processor ranks to create subgroups - INTEGER, ALLOCATABLE :: iwork(:,:), indx(:), isendcnts(:) ! see comments where used - INTEGER, ALLOCATABLE :: irecvcnts(:) ! see comments where used - INTEGER, ALLOCATABLE :: irecvcol(:) ! see comments where used - INTEGER, ALLOCATABLE :: ihavecolnum(:) ! Which columns do I have after first alltoall - INTEGER, ALLOCATABLE :: istartrecvbuf(:) ! Start of cols from proc i in recvbuf - INTEGER, ALLOCATABLE :: iinsendbuf(:) ! Now in sendbuf for proc i - INTEGER :: lensendbuf, lenrecvbuf ! Length of complete send- and recvbuf, first alltoall - INTEGER :: lenrecvbuf2,lensendbuf2 ! Length of complete send- and recvbuf, second alltoall - INTEGER :: jglob, jloc ! Global and local column index - INTEGER :: iglob, iloc ! Global and local row index - INTEGER :: hisjloc ! Local Achi column number of global column jglob sent - INTEGER :: numcols ! number of columns I get in first alltoall - INTEGER :: ierr1 ! For all the MPI calls - ! - ! Arrays to redistribute - ! - REAL,ALLOCATABLE :: sendbuffer_r(:) - REAL,ALLOCATABLE :: recvbuffer_r(:) - COMPLEX,ALLOCATABLE :: sendbuffer_c(:) - COMPLEX,ALLOCATABLE :: recvbuffer_c(:) - - LOGICAL:: l_real - l_real=PRESENT(achi_r) - ! - ! Now do the computation - ! - ierr=0 - np=nprow*npcol ! number of processors - mynumcols=m/np ! my number of columns in the column-cyclic distribution - IF (iam.LT.m-(m/np)*np) mynumcols=mynumcols+1 - myrow=iam/npcol ! my row number in the BLACS nprow*npcol grid - mycol=iam -(iam/npcol)*npcol ! my column number in the BLACS nprow*npcol grid - ! - ! Look, how many rows and columns Asca has locally - ! - myrowssca=(n-1)/(nb*nprow)*nb+ & - MIN(MAX(n-(n-1)/(nb*nprow)*nb*nprow-nb*myrow,0),nb) - IF (myrowssca.GT.lda) THEN - WRITE(6,*) 'Redist2: Wrong dimension of Asca, lda=',lda, & - ' myrowssca=',myrowssca - ENDIF - mycolssca=(m-1)/(nb*npcol)*nb+ & - MIN(MAX(m-(m-1)/(nb*npcol)*nb*npcol-nb*mycol,0),nb) - ! - ! Allocate arrays to build subgroups for rows and columns - ! and create subgroups consisting of processor rows and columns - ! - CALL MPI_COMM_GROUP(SUB_COMM,iworldgroup, Ierr) - IF (ierr.NE.0) THEN - WRITE(6,*) 'Mpi_comm_group Failed, ierr=',ierr - !CALL CPP_flush(6) - CALL juDFT_error('MPI_COMM_GROUP failed') - - ENDIF - ALLOCATE(iranks(1:MAX(npcol,nprow))) - ALLOCATE(irowgroup(0:nprow-1)) - ALLOCATE(icommrow(0:nprow-1)) - ALLOCATE(icolgroup(0:npcol-1)) - ALLOCATE(icommcol(0:npcol-1)) - DO k=0,nprow-1 ! Create groups of one processor row each - DO i=1,npcol - iranks(i)=k*npcol+i-1 - ENDDO - CALL MPI_GROUP_INCL(iworldgroup,npcol,iranks,irowgroup(k),& - ierr) - IF (ierr.NE.0) THEN - WRITE(6,*) 'MPI_GROUP_INCL(',k,') row failed, ierr=',ierr - !CALL CPP_flush(6) - CALL juDFT_error("MPI_GROUP_INCL failed") - ENDIF - CALL MPI_COMM_CREATE(SUB_COMM,irowgroup(k),& - icommrow(k),ierr) - IF (ierr.NE.0) THEN - WRITE(6,*) 'MPI_COMM_CREATE(',k,') row failed, ierr=',ierr - !CALL CPP_flush(6) - CALL juDFT_error("MPI_COMM_CREATE failed") - ENDIF - ENDDO - DO k=0,npcol-1 ! Create groups of one processor column each - DO i=1,nprow - iranks(i)=k+(i-1)*npcol - ENDDO - CALL MPI_GROUP_INCL(iworldgroup,nprow,iranks,icolgroup(k),& - ierr) - IF (ierr.NE.0) THEN - WRITE(6,*) 'MPI_GROUP_INCL(',k,') col failed, ierr=',ierr - !CALL CPP_flush(6) - CALL juDFT_error("MPI_GROUP_INCL failed") - ENDIF - CALL MPI_COMM_CREATE(SUB_COMM,icolgroup(k),& - icommcol(k),ierr) - IF (ierr.NE.0) THEN - WRITE(6,*) 'MPI_COMM_CREATE(',k,') col failed, ierr=',ierr - !CALL CPP_flush(6) - CALL juDFT_ERROR("MPI_COMM_CREATE failed") - ENDIF - ENDDO - DEALLOCATE(iranks) - ! - ! Now go through the local columns of Asca and look, to which - ! processor column they will belong - ! - ALLOCATE(iwork(1:mycolssca,0:MAX(npcol-1,nprow-1))) - ! iwork(i,jproc) will tell which local ScaLAPACK column is the - ! i-th column that will be sent to processor (myrow,modulo(jproc,npcol)) - ! in the first alltoallv - iwork=0 - ALLOCATE(indx(0:MAX(npcol-1,nprow-1))) - ! indx(modulo(jproc,npcol)) will indicate how many columns of processor column - ! modulo(jproc,npcol) are already found on this processor - indx=0 ! I have not found any column of any other processor yet - k=0 - DO jloc=1,mycolssca,nb - DO j=0,MIN(nb-1,mycolssca-jloc) - jglob=k*nb*npcol+mycol*nb+J+1 ! global column jglob, local jloc+j - jproc=MODULO(jglob-1,np) ! belonging to processor jproc - IF (jproc == iam) THEN - ! I will get this column on achi - ! There are (jglob-1)/np columns before this one on achi, and - ! It starts with element myrow*nb+1 and goes through it by steps - ! of nb*nprow - DO l=0,myrowssca/nb ! Go through column jglob, store to Achi - IF (l_real) THEN - DO i=1,MIN(nb,myrowssca-l*nb) - achi_r((jglob-1)/np*n+l*nprow*nb+myrow*nb+i)= & - asca_r(l*nb+i,jloc+j) - ENDDO - ELSE - DO i=1,MIN(nb,myrowssca-l*nb) - achi_c((jglob-1)/np*n+l*nprow*nb+myrow*nb+i)= & - asca_c(l*nb+i,jloc+j) - ENDDO - ENDIF - ENDDO - ELSE IF (MODULO(jproc,npcol) /= mycol) THEN - ! proc modulo(jglob-1,np) is not in my processor column, - ! but in column modulo(jproc,npcol) - ! I will send to him in first alltoallv - indx(MODULO(jproc,npcol))=indx(MODULO(jproc,npcol))+1 - iwork(indx(MODULO(jproc,npcol)),MODULO(jproc,npcol))= & - jloc+j - ! - ! I will send local column jloc+j to processor myrow,modulo(jproc,npcol) - ENDIF - ENDDO - k=k+1 ! dont forget, one local block is treated - ENDDO - ! Now I know what to send in first alltoallv: - ALLOCATE(isendcnts(0:MAX(npcol-1,nprow-1))) ! For alltoallv - isendcnts=0 - lensendbuf=0 ! sendbuffer is empty - DO k=0,npcol-1 - IF (k.NE.mycol) THEN - isendcnts(k)=indx(k)*myrowssca - lensendbuf=lensendbuf+isendcnts(k) - ENDIF - ENDDO - IF (l_real) THEN - ALLOCATE(sendbuffer_r(lensendbuf)) - ELSE - ALLOCATE(sendbuffer_c(lensendbuf)) - ENDIF - ! Put elements of asca to sendbuffer - ! - j=0 - DO k=0,npcol-1 - IF (k.NE.mycol) THEN - DO l=1,indx(k) - IF (l_real) THEN - DO i=1,myrowssca - sendbuffer_r(j+i)=asca_r(i,iwork(l,k)) - ENDDO - ELSE - DO i=1,myrowssca - sendbuffer_c(j+i)=asca_c(i,iwork(l,k)) - ENDDO - END IF - j=j+myrowssca - ENDDO - ENDIF - ENDDO - ! - ! Now look what I get in the first alltoallv - ! Go through the other processors in my row and see what they have for me - ! - ALLOCATE(irecvcnts(0:MAX(npcol-1,nprow-1))) - ALLOCATE(irecvcol(1:m)) - irecvcnts=0 - irecvcol=0 - lenrecvbuf=0 - numcols=0 ! Number of columns I will get in this alltoall - DO k=0,npcol-1 ! Now look what processors in my row send to me - IF (k.NE.mycol) THEN ! All except me are interesting - hiscolssca=(m-1)/(nb*npcol)*nb+ & - MIN(MAX(m-(m-1)/(nb*npcol)*nb*npcol-nb*k,0),nb) - ! His number of columns in ScaLPACK distribution - l=0 - DO jloc=1,hiscolssca,nb ! Blocks of nb columns - DO j=0,MIN(nb-1,hiscolssca-jloc) - jglob=l*nb*npcol+k*nb+J+1 ! global column jglob, local jloc+j - jproc=MODULO(jglob-1,np) ! Belonging to processor jproc for achi - IF ( MODULO(jproc,npcol).EQ.mycol) THEN - ! This proc is in my proc column - ! I will get this column in first alltoallv - numcols=numcols+1 ! I get one more column - irecvcol(numcols)=jglob ! it is global column - irecvcnts(k)=irecvcnts(k)+myrowssca ! more elements to rec - ENDIF - ENDDO - l=l+1 ! dont forget, one local block is treated - ENDDO - ENDIF ! If it is not me - lenrecvbuf=lenrecvbuf+irecvcnts(k) - ENDDO ! Loop over procs in my row - ! - ! Allocate receivebuffer - ! - IF (l_real) THEN - ALLOCATE(recvbuffer_r(lenrecvbuf)) - ELSE - ALLOCATE(recvbuffer_c(lenrecvbuf)) - END IF - ALLOCATE(iinsendbuf(0:MAX(npcol-1,nprow-1))) - iinsendbuf(0)=0 - ALLOCATE(istartrecvbuf(0:MAX(npcol-1,nprow-1))) - istartrecvbuf(0)=0 - DO l=0,npcol-2 - IF (l.EQ.mycol) THEN - iinsendbuf(l+1)=iinsendbuf(l) - istartrecvbuf(l+1)=istartrecvbuf(l) - ELSE - iinsendbuf(l+1)=iinsendbuf(l)+isendcnts(l) - istartrecvbuf(l+1)=istartrecvbuf(l)+irecvcnts(l) - ENDIF - ENDDO - IF(l_real) THEN - CALL MPI_ALLTOALLV(sendbuffer_r,isendcnts,iinsendbuf,& - CPP_MPI_TYP_REAL,recvbuffer_r,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_REAL,icommrow(myrow),ierr) - ELSE - CALL MPI_ALLTOALLV(sendbuffer_c,isendcnts,iinsendbuf,& - CPP_MPI_TYP_COMPLEX,recvbuffer_c,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_COMPLEX,icommrow(myrow),ierr) - ENDIF - ! - ! After first alltoallv deallocate sendbuffer (not necessary) - ! And set indexbuffers of sendbuffer to zero - ! - ! deallocate(sendbuffer) - lensendbuf2=0 - iinsendbuf(0)=0 - isendcnts=0 - ! - ! Now look how many columns of Achi the processors in my processor - ! column have. I should have myrossca elements of all of these rows - ! either in my Asca or in my receivebuffer - ! - indx=m/np ! Each processor has at least m/np columns - DO k=0,nprow-1 ! Procs in my column - ! global processor number is k*npcol+mycol, first ones possibly one more column - IF (k*npcol+mycol.LT.m-(m/np*np)) indx(k)=indx(k)+1 - IF (k.NE.myrow) THEN - isendcnts(k)=indx(k)*myrowssca - lensendbuf2=lensendbuf2+isendcnts(k) - IF (k.LT.nprow-1) iinsendbuf(k+1)=iinsendbuf(k)+isendcnts(k) - ELSE - isendcnts(myrow)=0 - IF (k.LT.nprow-1) iinsendbuf(k+1)=iinsendbuf(k) - ENDIF - ENDDO ! End k-loop over procs in my column - IF (lensendbuf2 > lensendbuf) THEN - IF (l_real) THEN - DEALLOCATE(sendbuffer_r) - ALLOCATE(sendbuffer_r(lensendbuf2)) - ELSE - DEALLOCATE(sendbuffer_c) - ALLOCATE(sendbuffer_c(lensendbuf2)) - END IF - - ENDIF - ! - ! Now look through Asca which columns will be sent now - ! And put them to the correct place in sendbuf - ! - k=0 - DO jloc=1,mycolssca,nb - DO j=0,MIN(nb-1,mycolssca-jloc) - jglob=k*nb*npcol+mycol*nb+J+1 ! global column jglob, local jloc+j - jproc=MODULO(jglob-1,np) ! belonging to processor jproc - ! if (jproc.eq.iam) then - ! - ! Everything done before first alltoall, columns already on achi - ! - ! else if (modulo(jproc,npcol).ne.mycol) then - ! - ! proc modulo(jglob-1,np) is not in my processor column, - ! I sent it to him in first alltoallv - ! - IF (MODULO(jproc,npcol) == mycol.AND.jproc /= iam) THEN - ! - ! proc modulo(jglob-1,np) is in my column, I will send to him - ! now in second alltoallv - ! - jprocrow=MODULO(jglob-1,np)/npcol ! and in this row - ! Now look which local column number for this proc column - ! jglob will be and put it to the correct place in sendbuffer - ! This place is iinsendbuf(k)+ - ! number of local columns of proc jprocrow,mycol before column jglob - ! * myrowssca + i,i=1,myrowssca - hisjloc=(jglob-1)/np+1 ! his local Achi column number - ! - ! Store to sendbuffer - ! - IF (l_real) THEN - DO i=1,myrowssca - sendbuffer_r(iinsendbuf(jprocrow)+(hisjloc-1)*myrowssca+i) =asca_r(i,jloc+j) - ENDDO - ELSE - DO i=1,myrowssca - sendbuffer_c(iinsendbuf(jprocrow)+(hisjloc-1)*myrowssca+i) =asca_c(i,jloc+j) - ENDDO - ENDIF - ENDIF - ENDDO - k=k+1 ! dont forget, one local block is treated - ENDDO - ! - ! and through receivebuffer which columns are now there and - ! which will be sent on - ! - DO k=1,numcols !number of columns I got - jglob=irecvcol(k) ! The global column number of column k in the recvbuf - jproc=MODULO(jglob-1,np) ! The processor it will belong to is jproc - IF (jproc.EQ.iam) THEN - ! - ! It is my column, have to store it to achi - ! - jloc=(jglob-1)/np+1 ! Local Achi column number - l=0 ! Number of blocks already treated - DO i=1,myrowssca,nb ! Go through row blocks - IF (l_real) THEN - DO kk=0,MIN(nb-1,myrowssca-i) - achi_r((jloc-1)*n+myrow*nb+l*nb*nprow+kk+1) =recvbuffer_r((k-1)*myrowssca+i+kk) - ENDDO - ELSE - DO kk=0,MIN(nb-1,myrowssca-i) - achi_c((jloc-1)*n+myrow*nb+l*nb*nprow+kk+1) =recvbuffer_c((k-1)*myrowssca+i+kk) - ENDDO - ENDIF - l=l+1 ! Next row block - ENDDO - ELSE ! It is not my column, must be sent to jprocrow - jprocrow=MODULO(jglob-1,np)/npcol ! jproc's row number ' - jloc=(jglob-1)/np+1 ! Local Achi column number for jproc - ! - ! Store it to sendbuffer - ! - IF (l_real) THEN - DO i=1,myrowssca - sendbuffer_r(iinsendbuf(jprocrow)+(jloc-1)*myrowssca+i) & - =recvbuffer_r((k-1)*myrowssca+i) - ENDDO - ELSE - DO i=1,myrowssca - sendbuffer_c(iinsendbuf(jprocrow)+(jloc-1)*myrowssca+i) & - =recvbuffer_c((k-1)*myrowssca+i) - ENDDO - ENDIF - ENDIF - ENDDO - ! - ! Now look how large the next receivebuffer has to be - ! - lenrecvbuf2=0 - irecvcnts=0 - istartrecvbuf=0 - DO k=0,nprow-1 ! Procs in my column - IF (k.NE.myrow) THEN ! I get from this proc - hisrowssca= (n-1)/(nb*nprow)*nb+ & - MIN(MAX(n-(n-1)/(nb*nprow)*nb*nprow-nb*k,0),nb) ! Proc k,mycol has - ! hisrowssca rows in ScaLAPACK distribution - ! It sends indx(myrow) columns of that length because it has elements of - ! all of my columns - irecvcnts(k)=indx(myrow)*hisrowssca ! I'll get irecvcnts(k) elements from ' - ! ( k, mycol) - lenrecvbuf2=lenrecvbuf2+irecvcnts(k) - IF (k.LT.nprow-1) istartrecvbuf(k+1)= & - istartrecvbuf(k)+irecvcnts(k) - ELSE ! Nothing to receive, it is me - irecvcnts(myrow)=0 - IF (k.LT.nprow-1) istartrecvbuf(k+1)=istartrecvbuf(k) - ENDIF - ENDDO ! End k-loop over procs in my column - IF (lenrecvbuf2 > lenrecvbuf) THEN - IF (l_real) THEN - DEALLOCATE(recvbuffer_r) - ALLOCATE(recvbuffer_r(lenrecvbuf2)) - ELSE - DEALLOCATE(recvbuffer_c) - ALLOCATE(recvbuffer_c(lenrecvbuf2)) - ENDIF - ENDIF - ! - IF (l_real) THEN - CALL MPI_ALLTOALLV(sendbuffer_r,isendcnts,iinsendbuf,& - CPP_MPI_TYP_REAL,recvbuffer_r,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_REAL,icommcol(mycol),ierr) - ELSE - CALL MPI_ALLTOALLV(sendbuffer_c,isendcnts,iinsendbuf,& - CPP_MPI_TYP_COMPLEX,recvbuffer_c,irecvcnts,istartrecvbuf,& - CPP_MPI_TYP_COMPLEX,icommcol(mycol),ierr) - ENDIF - ! - ! After second alltoallv deallocate sendbuffer - ! And indexbuffers - ! - DEALLOCATE(iwork) - DEALLOCATE(isendcnts) - DEALLOCATE(iinsendbuf) - ! - ! Store the received elements to Achi - ! - DO k=0,nprow-1 ! All except me have sent something - IF (k.NE.myrow) THEN - hisrowssca= (n-1)/(nb*nprow)*nb+ & - MIN(MAX(n-(n-1)/(nb*nprow)*nb*nprow-nb*k,0),nb) - DO j=1,indx(myrow) ! Go through my columns, j is local Achi col num - l=0 ! Number of blocks of this column already treated - DO i=1,hisrowssca,nb ! Go through his rowblocks - IF (l_real) THEN - DO kk=0,MIN(nb-1,hisrowssca-i) - achi_r((j-1)*n+k*nb+l*nb*nprow+kk+1) =recvbuffer_r(istartrecvbuf(k)+(j-1)*hisrowssca+i+kk) - ENDDO - ELSE - DO kk=0,MIN(nb-1,hisrowssca-i) - achi_c((j-1)*n+k*nb+l*nb*nprow+kk+1) =recvbuffer_c(istartrecvbuf(k)+(j-1)*hisrowssca+i+kk) - ENDDO - ENDIF - l=l+1 ! Next row block - ENDDO - ENDDO - ENDIF - ENDDO - DEALLOCATE(indx) - DEALLOCATE(irecvcnts) - DEALLOCATE(istartrecvbuf) - ! - ! Free the row and column communicators - ! - DO k=0,nprow-1 - IF (icommrow(k).NE.MPI_COMM_NULL) THEN - CALL MPI_COMM_FREE(icommrow(k),ierr) - ENDIF - CALL MPI_GROUP_FREE(irowgroup(k),ierr) - ENDDO - DO k=0,npcol-1 - IF (icommcol(k).NE.MPI_COMM_NULL) THEN - CALL MPI_COMM_FREE(icommcol(k),ierr) - ENDIF - CALL MPI_GROUP_FREE(icolgroup(k),ierr) - ENDDO - ! - ! Deallocate the arrays to store the indices - ! - DEALLOCATE(irecvcol) - DEALLOCATE(irowgroup) - DEALLOCATE(icommrow) - DEALLOCATE(icolgroup) - DEALLOCATE(icommcol) - RETURN -END SUBROUTINE subredist2 -END diff --git a/diagonalization/zsymsecloc.F90 b/diagonalization/zsymsecloc.F90 deleted file mode 100644 index cdb59909205a61db071244b84e7c534a261ed5f5..0000000000000000000000000000000000000000 --- a/diagonalization/zsymsecloc.F90 +++ /dev/null @@ -1,41 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -! NOTE: this contains only the interface, the actual code is included by the -! proprocessor from the file zsymsecloc_cpp.F90 -! - -MODULE m_zsymsecloc - use m_juDFT -!******************************************************* -! Solve the generalized secular equation. -! For film-systems exhibiting -! z-reflexion symmetry, the basis is transformed to -! even and odd functions and the even-even and odd-odd -! blocks are diagonalized separately. -! If local orbitals are present in a film with z-reflection, -! locrectify is used to construct linear combinations of -! the local orbitals that are eigenfunctions of the z- -! reflexion operation. -! Frank Freimuth, January 2006 -!******************************************************* - INTERFACE zsymsecloc - MODULE procedure zsymsecloc_r,zsymsecloc_c - END INTERFACE zsymsecloc -CONTAINS - SUBROUTINE zsymsecloc_r(jsp,input,lapw,bkpt,atoms, kveclo, sym,l_zref,cell, dimension,matsize, nsize, jij,matind,nred,eig,ne, a,b, z) - -#define CPP_REALDATA -#include "zsymsecloc_cpp.F90" - END SUBROUTINE zsymsecloc_r - - SUBROUTINE zsymsecloc_c(jsp,input,lapw,bkpt,atoms, kveclo, sym,l_zref,cell, dimension,matsize, nsize, jij,matind,nred,eig,ne, a,b, z) - -#undef CPP_REALDATA -#include "zsymsecloc_cpp.F90" - END SUBROUTINE zsymsecloc_c - -END MODULE m_zsymsecloc diff --git a/diagonalization/zsymsecloc_cpp.F90 b/diagonalization/zsymsecloc_cpp.F90 deleted file mode 100644 index 860e8f68ca1cef3fe8747ffe35fe6569dd186507..0000000000000000000000000000000000000000 --- a/diagonalization/zsymsecloc_cpp.F90 +++ /dev/null @@ -1,518 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -! This file is included in zsymsecloc.F90 twice with CPP_REALDATA set and unset -! to create both code for real and complex - -#include"cpp_double.h" - - USE m_locrectify - USE m_geneigprobl - USE m_types - IMPLICIT NONE - TYPE(t_dimension),INTENT(IN) :: dimension - TYPE(t_input),INTENT(IN) :: input - TYPE(t_jij),INTENT(IN) :: jij - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_atoms),INTENT(IN) :: atoms - - TYPE(t_lapw),INTENT(IN) :: lapw - real,intent(in) ::bkpt(3) - integer,intent(in)::kveclo(atoms%nlotot) - logical,intent(in) :: l_zref - - ! .. - ! .. Scalar Arguments .. - INTEGER, INTENT (IN) :: matsize,jsp - INTEGER, INTENT (INOUT) :: nsize - INTEGER, INTENT (OUT) :: ne - INTEGER, INTENT (IN) :: nred - ! .. - ! .. Array Arguments .. - INTEGER, INTENT (IN) :: matind(dimension%nbasfcn,2) - REAL, INTENT (OUT) :: eig(dimension%neigd) - -#ifdef CPP_REALDATA - REAL, ALLOCATABLE, INTENT (INOUT) :: a(:),b(:) - REAL, ALLOCATABLE, INTENT (INOUT) :: z(:,:) -#else - COMPLEX, ALLOCATABLE, INTENT (INOUT)::a(:),b(:) - COMPLEX, ALLOCATABLE, INTENT (INOUT) :: z(:,:) -#endif - - -#ifdef CPP_REALDATA - real locrec(atoms%nlotot,atoms%nlotot) - real,allocatable::atemp(:,:) - real,allocatable::btemp(:,:) - real,allocatable::aa(:) - real,allocatable::bb(:) - real,allocatable::z1(:,:) - real,allocatable::z2(:,:) - real,allocatable::ztemp(:,:) - real recsqrtwo,dobzero -#else - complex locrec(atoms%nlotot,atoms%nlotot) - complex,allocatable:: atemp(:,:) - complex,allocatable:: btemp(:,:) - complex,allocatable::aa(:) - complex,allocatable::bb(:) - complex,allocatable::z1(:,:) - complex,allocatable::z2(:,:) - complex,allocatable::ztemp(:,:) - complex recsqrtwo,dobzero -#endif - logical l_verbose - logical kindlocrec(atoms%nlotot) - real,allocatable::etemp1(:),etemp2(:) - integer iind,jind,ii,info,iu,pos1,pos2,pos3,i2,j2 - integer ne1,ne2,i1,j1,i,jg,j - logical,allocatable:: evensort(:) - integer evenlocs,oddlocs - integer evenlocindex(atoms%nlotot) - integer oddlocindex(atoms%nlotot) - - - ! print*,"in zsymsecloc" - deallocate(z) - - !****************************************** - ! l_zref=.false. => simply call eigensolver - !****************************************** - if(.not.l_zref)then -#ifdef CPP_REALDATA - call geneigprobl(dimension%nbasfcn, nsize,dimension%neigd,jij%l_j,eig,ne,a,b,z) -#else - call geneigprobl(dimension%nbasfcn, nsize,dimension%neigd,jij%l_j,eig,ne,a_c=a,b_c=b,z_c=z) -#endif - allocate(a(dimension%nbasfcn*(dimension%nbasfcn+1)/2)) - allocate(b(dimension%nbasfcn*(dimension%nbasfcn+1)/2)) - return - !****************************************** - ! l_zref=.true. => blockdiagonalize - ! hamiltonian and overlap-matrix by - ! transforming to even/odd basisfunctions - !****************************************** - else - inquire(file='verbose',exist=l_verbose) - if(.not.l_verbose)inquire(file='veryverbose',exist=l_verbose) - iu=dimension%neigd/2 -#ifdef CPP_REALDATA - recsqrtwo=1.0/sqrt(2.0) - dobzero=0.0 -#else - recsqrtwo=cmplx(1.0/sqrt(2.0),0.0) - dobzero=cmplx(0.0,0.0) -#endif - - - if(atoms%nlotot.gt.0)then - !********************************************* - ! Find even/odd linear combinations of locs - !********************************************* - if(l_verbose)then - print*,"find lincos of locs that are eigenfunctions of zreflection" - print*,"apws=",lapw%nv(jsp) - print*,"atoms%nlotot=",atoms%nlotot - print*,"basis-size=",nsize - endif - IF(nsize/=(lapw%nv(jsp)+atoms%nlotot)) & - CALL juDFT_error("nsize /= lapw%nv + atoms%nlotot" ,calledby ="zsymsecloc") -#ifdef CPP_REALDATA - call locrectify(jsp,input,lapw,bkpt,atoms, kveclo, sym,cell,& - kindlocrec,evenlocs,oddlocs, evenlocindex,oddlocindex,locrec) -#else - call locrectify(jsp,input,lapw,bkpt,atoms, kveclo, sym,cell,& - kindlocrec,evenlocs,oddlocs, evenlocindex,oddlocindex,locrec_c=locrec) -#endif - - - !********************************************* - ! Perform basis-transformation of Hamiltonian - ! and Overlap matrix. The new basis is the basis - ! of even / odd (with respect to z-reflection) - ! local orbitals. - !********************************************* - allocate(atemp(lapw%nv(jsp),atoms%nlotot)) - allocate(btemp(atoms%nlotot,lapw%nv(jsp))) - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - atemp(1:lapw%nv(jsp),iind)=a(pos1+1:pos1+lapw%nv(jsp)) - enddo -#ifdef CPP_REALDATA - call CPP_BLAS_sgemm('T','T',atoms%nlotot,lapw%nv(jsp),atoms%nlotot,real(1.0),& - locrec,atoms%nlotot,atemp,lapw%nv(jsp),real(0.0),btemp,atoms%nlotot) -#else - call CPP_BLAS_cgemm('C','T',atoms%nlotot,lapw%nv(jsp),atoms%nlotot,cmplx(1.0,0.0),& - locrec,atoms%nlotot,atemp,lapw%nv(jsp),cmplx(0.0,0.0),btemp,atoms%nlotot) -#endif - atemp(:,:) = transpose(btemp) - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - a(pos1+1:pos1+lapw%nv(jsp))=atemp(1:lapw%nv(jsp),iind) - enddo - - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - atemp(1:lapw%nv(jsp),iind)=b(pos1+1:pos1+lapw%nv(jsp)) - enddo -#ifdef CPP_REALDATA - call CPP_BLAS_sgemm('T','T',atoms%nlotot,lapw%nv(jsp),atoms%nlotot,real(1.0),& - locrec,atoms%nlotot,atemp,lapw%nv(jsp),real(0.0),btemp,atoms%nlotot) -#else - call CPP_BLAS_cgemm('C','T',atoms%nlotot,lapw%nv(jsp),atoms%nlotot,cmplx(1.0,0.0),& - locrec,atoms%nlotot,atemp,lapw%nv(jsp),cmplx(0.0,0.0),btemp,atoms%nlotot) -#endif - atemp(:,:) = transpose(btemp) - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - b(pos1+1:pos1+lapw%nv(jsp))=atemp(1:lapw%nv(jsp),iind) - enddo - deallocate(atemp) - deallocate(btemp) - - allocate(atemp(atoms%nlotot,atoms%nlotot)) - allocate(btemp(atoms%nlotot,atoms%nlotot)) - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - do jind=1,iind-1 - atemp(iind,jind)=a(pos1+lapw%nv(jsp)+jind) -#ifdef CPP_REALDATA - atemp(jind,iind)=a(pos1+lapw%nv(jsp)+jind) -#else - atemp(jind,iind)=conjg(a(pos1+lapw%nv(jsp)+jind)) -#endif - enddo - atemp(iind,iind)=a(pos1+lapw%nv(jsp)+iind) - enddo -#ifdef CPP_REALDATA - call CPP_BLAS_sgemm('T','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,real(1.0),& - locrec,atoms%nlotot,atemp,atoms%nlotot,real(0.0),btemp,atoms%nlotot) - call CPP_BLAS_sgemm('N','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,real(1.0),& - btemp,atoms%nlotot,locrec,atoms%nlotot,real(0.0),atemp,atoms%nlotot) -#else - call CPP_BLAS_cgemm('C','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,cmplx(1.0,0.0),& - locrec,atoms%nlotot,atemp,atoms%nlotot,cmplx(0.0,0.0),btemp,atoms%nlotot) - call CPP_BLAS_cgemm('N','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,cmplx(1.0,0.0),& - btemp,atoms%nlotot,locrec,atoms%nlotot,cmplx(0.0,0.0),atemp,atoms%nlotot) -#endif - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - do jind=1,iind - a(pos1+lapw%nv(jsp)+jind)=atemp(iind,jind) - enddo - enddo - - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - do jind=1,iind-1 - atemp(iind,jind)=b(pos1+lapw%nv(jsp)+jind) -#ifdef CPP_REALDATA - atemp(jind,iind)=b(pos1+lapw%nv(jsp)+jind) -#else - atemp(jind,iind)=conjg(b(pos1+lapw%nv(jsp)+jind)) -#endif - enddo - atemp(iind,iind)=b(pos1+lapw%nv(jsp)+iind) - enddo -#ifdef CPP_REALDATA - call CPP_BLAS_sgemm('T','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,real(1.0),& - locrec,atoms%nlotot,atemp,atoms%nlotot,real(0.0),btemp,atoms%nlotot) - call CPP_BLAS_sgemm('N','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,real(1.0),& - btemp,atoms%nlotot,locrec,atoms%nlotot,real(0.0),atemp,atoms%nlotot) -#else - call CPP_BLAS_cgemm('C','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,cmplx(1.0,0.0),& - locrec,atoms%nlotot,atemp,atoms%nlotot,cmplx(0.0,0.0),btemp,atoms%nlotot) - call CPP_BLAS_cgemm('N','N',atoms%nlotot,atoms%nlotot,atoms%nlotot,cmplx(1.0,0.0),& - btemp,atoms%nlotot,locrec,atoms%nlotot,cmplx(0.0,0.0),atemp,atoms%nlotot) -#endif - do iind=1,atoms%nlotot - pos1=((lapw%nv(jsp)+iind)*(lapw%nv(jsp)+iind-1))/2 - do jind=1,iind - b(pos1+lapw%nv(jsp)+jind)=atemp(iind,jind) - enddo - enddo - - deallocate(atemp) - deallocate(btemp) - - else - evenlocs=0 - oddlocs=0 - endif ! atoms%nlotot.gt.0 - - !********************************************* - ! Test matind. - !********************************************* - ii=0 - jind=0 - pos1=0 - do iind=1,nred - IF(matind(iind,1) eigenvalue ii belongs to even spectrum - !evensort(ii)=.false. => eigenvalue ii belongs to odd spectrum - do ii=1,ne - if(etemp1(iind).lt.etemp2(jind)) then - evensort(ii)=.true. - iind=iind+1 - else - evensort(ii)=.false. - jind=jind+1 - endif - enddo - iind=1 !position in the even-arrays - jind=1 !position in the oneD%odd-arrays - do ii=1,ne - if(evensort(ii))then - eig(ii)=etemp1(iind) - z(1:lapw%nv(jsp)+atoms%nlotot,ii)=dobzero - !Recover the eigenvectors of the original problem for the even block - do i=1,nred - z(matind(i,1),ii)=z1(i,iind)*recsqrtwo - z(matind(i,2),ii)=z(matind(i,2),ii)+z1(i,iind)*recsqrtwo - enddo !i - if(atoms%nlotot.gt.0) z(lapw%nv(jsp)+1:lapw%nv(jsp)+atoms%nlotot,ii)=dobzero - do pos1=1,evenlocs - do i=1,atoms%nlotot -#ifdef CPP_REALDATA - z(lapw%nv(jsp)+i,ii)=z(lapw%nv(jsp)+i,ii)+ z1(nred+pos1,iind)*locrec(i,evenlocindex(pos1)) -#else - z(lapw%nv(jsp)+i,ii)=z(lapw%nv(jsp)+i,ii)+ z1(nred+pos1,iind)*conjg(locrec(i,evenlocindex(pos1))) -#endif - enddo - enddo - iind=iind+1 - else - !Recover the eigenvectors of the original problem for the odd block - eig(ii)=etemp2(jind) - j1=0 - do i=1,nred - if(matind(i,1).ne.matind(i,2))then - j1=j1+1 - z(matind(i,1),ii)=z2(j1,jind)*recsqrtwo - z(matind(i,2),ii)=-z2(j1,jind)*recsqrtwo - else - z(matind(i,1),ii)=dobzero - endif - enddo !i - if(atoms%nlotot.gt.0) z(lapw%nv(jsp)+1:lapw%nv(jsp)+atoms%nlotot,ii)=dobzero - do pos1=1,oddlocs - do i=1,atoms%nlotot -#ifdef CPP_REALDATA - z(lapw%nv(jsp)+i,ii)=z(lapw%nv(jsp)+i,ii)+ z2(lapw%nv(jsp)-nred+pos1,jind)*locrec(i,oddlocindex(pos1)) -#else - z(lapw%nv(jsp)+i,ii)=z(lapw%nv(jsp)+i,ii)+ z2(lapw%nv(jsp)-nred+pos1,jind)*conjg(locrec(i,oddlocindex(pos1))) -#endif - enddo - enddo - jind=jind+1 - endif !evensort - enddo !ii - - allocate(a(dimension%nbasfcn*(dimension%nbasfcn+1)/2)) - allocate(b(dimension%nbasfcn*(dimension%nbasfcn+1)/2)) - endif !l_zref - - deallocate ( z1,z2,etemp1,etemp2,evensort ) diff --git a/eigen/soc_init.f90 b/eigen/soc_init.f90 deleted file mode 100644 index b95e25432496378616b8d64a0a259c131bae867d..0000000000000000000000000000000000000000 --- a/eigen/soc_init.f90 +++ /dev/null @@ -1,116 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -MODULE m_socinit - USE m_juDFT - IMPLICIT NONE - - -CONTAINS - !>Initialization of SOC matrix elements used in first variation SOC - !! - SUBROUTINE socinit(mpi,atoms,sphhar,enpara,input,vr,noco,& !in - usdus,rsoc) !out - !Initialized the radial-spin-orbit elements in rsoc - !needed for first variation SOC - USE m_soinit - USE m_types - IMPLICIT NONE - TYPE(t_mpi),INTENT(IN) :: mpi - TYPE(t_input),INTENT(IN) :: input - TYPE(t_noco),INTENT(IN) :: noco - TYPE(t_sphhar),INTENT(IN) :: sphhar - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_enpara),INTENT(IN) :: enpara - ! .. - ! .. Scalar Arguments .. - ! .. - ! .. Array Arguments .. - REAL, INTENT (IN) :: vr(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd) - TYPE(t_usdus),INTENT(INOUT):: usdus - TYPE(t_rsoc),INTENT(OUT) :: rsoc - - ! .. - ! .. Local Scalars .. - INTEGER l,n - ! for Spin-orbit... - LOGICAL :: l_test - LOGICAL, SAVE :: first_k = .TRUE. - - - CHARACTER*3 :: chntype - - - ALLOCATE(rsoc%rsopp(atoms%ntype,atoms%lmaxd,2,2),rsoc%rsoppd (atoms%ntype,atoms%lmaxd,2,2)) - ALLOCATE(rsoc%rsopdp (atoms%ntype,atoms%lmaxd,2,2),rsoc%rsopdpd(atoms%ntype,atoms%lmaxd,2,2)) - ALLOCATE(rsoc%rsoplop (atoms%ntype,atoms%nlod,2,2),rsoc%rsoplopd(atoms%ntype,atoms%nlod,2,2)) - ALLOCATE(rsoc%rsopdplo(atoms%ntype,atoms%nlod,2,2),rsoc%rsopplo (atoms%ntype,atoms%nlod,2,2)) - ALLOCATE(rsoc%rsoploplop(atoms%ntype,atoms%nlod,atoms%nlod,2,2)) - - - CALL soinit(atoms,input,enpara,vr,noco%soc_opt(atoms%ntype+2),& - rsoc%rsopp,rsoc%rsoppd,rsoc%rsopdp,rsoc%rsopdpd,usdus,& - rsoc%rsoplop,rsoc%rsoplopd,rsoc%rsopdplo,rsoc%rsopplo,rsoc%rsoploplop) - INQUIRE(file="socscale",exist=l_test) - IF (l_test) THEN - OPEN(99,file="socscale") - READ(99,*) n - CLOSE(99) - WRITE(*,*) "SOC scaled by ",n,"%" - rsoc%rsopp(:,:,:,:) = n/100.* rsoc%rsopp - rsoc%rsopdp(:,:,:,:) = n/100.*rsoc%rsopdp - rsoc%rsoppd(:,:,:,:) = n/100.*rsoc%rsoppd - rsoc%rsopdpd(:,:,:,:) = n/100.*rsoc%rsopdpd - rsoc%rsoplop(:,:,:,:) = n/100.*rsoc%rsoplop - rsoc%rsoplopd(:,:,:,:) = n/100.*rsoc%rsoplopd - rsoc%rsopdplo(:,:,:,:) = n/100.*rsoc%rsopdplo - rsoc%rsopplo(:,:,:,:) = n/100.* rsoc%rsopplo - rsoc%rsoploplop(:,:,:,:,:) = n/100.*rsoc%rsoploplop - ENDIF - IF (noco%soc_opt(atoms%ntype+1)) THEN - DO n= 1,atoms%ntype - IF (.NOT. noco%soc_opt(n)) THEN - rsoc%rsopp(n,:,:,:) = 0.0 - rsoc%rsopdp(n,:,:,:) = 0.0 - rsoc%rsoppd(n,:,:,:) = 0.0 - rsoc%rsopdpd(n,:,:,:) = 0.0 - rsoc%rsoplop(n,:,:,:) = 0.0 - rsoc%rsoplopd(n,:,:,:) = 0.0 - rsoc%rsopdplo(n,:,:,:) = 0.0 - rsoc%rsopplo(n,:,:,:) = 0.0 - rsoc%rsoploplop(n,:,:,:,:) = 0.0 - ENDIF - ENDDO - ENDIF - - IF ((first_k).AND.(mpi%irank.EQ.0)) THEN - DO n = 1,atoms%ntype - WRITE (6,FMT=8000) - WRITE (6,FMT=9000) - WRITE (6,FMT=8001) (2*rsoc%rsopp(n,l,1,1),l=1,3) - WRITE (6,FMT=8001) (2*rsoc%rsopp(n,l,2,2),l=1,3) - WRITE (6,FMT=8001) (2*rsoc%rsopp(n,l,2,1),l=1,3) - ENDDO - IF (noco%soc_opt(atoms%ntype+1)) THEN - WRITE (chntype,'(i3)') atoms%ntype - WRITE (6,fmt='(A,2x,'//chntype//'l1)') 'SOC contribution of certain atom types included in Hamiltonian:',& - (noco%soc_opt(n),n=1,atoms%ntype) - ELSE - WRITE(6,fmt='(A,1x,A)') 'SOC contribution of all atom types included in Hamiltonian.' - ENDIF - IF (noco%soc_opt(atoms%ntype+2)) THEN - WRITE(6,fmt='(A)') 'SOC Hamiltonian is constructed by neglecting B_xc.' - ENDIF - first_k=.FALSE. - ENDIF -8000 FORMAT (' spin - orbit parameter HR ') -8001 FORMAT (8f8.4) -9000 FORMAT (5x,' p ',5x,' d ', 5x, ' f ') - - - RETURN - END SUBROUTINE socinit -END MODULE m_socinit diff --git a/eigen_soc/soinit.f90 b/eigen_soc/soinit.f90 deleted file mode 100644 index 00d2394af9d11ad06aac1436795f3420a2602f5a..0000000000000000000000000000000000000000 --- a/eigen_soc/soinit.f90 +++ /dev/null @@ -1,59 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- -MODULE m_soinit - ! - !********************************************************************** - ! generates radial spin-orbit matrix elements:sorad - !********************************************************************** - ! -CONTAINS - SUBROUTINE soinit(atoms,input,enpara, vr,spav,rsoc,usdus) - - USE m_sorad - USE m_types - IMPLICIT NONE - TYPE(t_enpara),INTENT(IN) :: enpara - TYPE(t_input),INTENT(IN) :: input - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_usdus),INTENT(INOUT):: usdus - TYPE(t_rsoc),INTENT(INOUT) :: rsoc - ! - ! .. Scalar Arguments .. - ! .. - ! .. Scalar Arguments .. - LOGICAL, INTENT (IN) :: spav ! if T, spin-averaged pot is used - ! .. - ! .. Array Arguments .. - REAL, INTENT (IN) :: vr(:,0:,:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) - ! .. - ! .. Local Scalars .. - INTEGER i,jspin,n - ! .. - ! .. Local Arrays .. - REAL vr0(atoms%jmtd,size(vr,4)) - ! .. - rsoc%rsopp =0.0 - rsoppd =0.0 - rsopdp =0.0 - rsopdpd=0.0 - rsoplop =0.0 - rsoplopd=0.0 - rsopdplo=0.0 - rsopplo=0.0 - rsoploplop=0.0 - DO n = 1,atoms%ntype - vr0=0.0 - vr0(:atoms%jri(n),:) = vr(:atoms%jri(n),0,n,:) - - CALL sorad(& - atoms,input,n,vr0,enpara,spav,& - rsopp,rsopdpd,rsoppd,rsopdp,usdus,& - rsoplop,rsoplopd,rsopdplo,rsopplo,rsoploplop) - - END DO ! end-do-loop : atoms%ntype - - END SUBROUTINE soinit -END MODULE m_soinit diff --git a/global/enpara.f90 b/global/enpara.f90 deleted file mode 100644 index be6629f2e51f83597957a88c4678e09671fb1638..0000000000000000000000000000000000000000 --- a/global/enpara.f90 +++ /dev/null @@ -1,366 +0,0 @@ - MODULE m_enpara - use m_juDFT - USE m_constants -! ************************************************************* -! Module containing three subroutines -! r_enpara: read enpara file -! w_enpara: write enpara file -! mix_enpara: calculate new energy parameters -! ************************************************************* - CONTAINS - SUBROUTINE w_enpara( & - & atoms,jspin,film,& - & enpara,& - & id) -! -! write enpara-file -! - USE m_types - IMPLICIT NONE - - INTEGER, INTENT (IN) :: jspin,id - LOGICAL,INTENT(IN) :: film - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_enpara),INTENT(IN) :: enpara - - INTEGER n,l,lo - LOGICAL l_opened - - INQUIRE(unit=40,OPENED=l_opened) - if (.not.l_opened) return - - - WRITE (40,FMT=8035) jspin,enpara%enmix(jspin) - WRITE (40,FMT=8036) - 8035 FORMAT (5x,'energy parameters for spin ',i1,' mix=',f10.6) - 8036 FORMAT (t6,'atom',t15,'s',t24,'p',t33,'d',t42,'f') - DO n = 1,atoms%ntype - WRITE (6,FMT=8040) n, (enpara%el0(l,n,jspin),l=0,3),& - & (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin) - WRITE (id,FMT=8040) n, (enpara%el0(l,n,jspin),l=0,3),& - & (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin) - WRITE (40,FMT=8040) n, (enpara%el0(l,n,jspin),l=0,3),& - & (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin) -!---> energy parameters for the local orbitals - IF (atoms%nlo(n).GE.1) THEN - WRITE (6,FMT=8039) (enpara%ello0(lo,n,jspin),lo=1,atoms%nlo(n)) - WRITE (6,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n)) - WRITE (id,FMT=8039) (enpara%ello0(lo,n,jspin),lo=1,atoms%nlo(n)) - WRITE (id,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n)) - WRITE (40,FMT=8039) (enpara%ello0(lo,n,jspin),lo=1,atoms%nlo(n)) - WRITE (40,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n)) - END IF - - ENDDO - 8038 FORMAT (' --> change ',60(l1,8x)) - 8039 FORMAT (' --> lo ',60f9.5) - 8040 FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3) - - IF (film) THEN - WRITE (40,FMT=8050) enpara%evac0(1,jspin),enpara%lchg_v(1,jspin),enpara%evac0(2,jspin) - WRITE (6,FMT=8050) enpara%evac0(1,jspin),enpara%lchg_v(1,jspin),enpara%evac0(2,jspin) - WRITE (id,FMT=8050) enpara%evac0(1,jspin),enpara%lchg_v(1,jspin),enpara%evac0(2,jspin) - 8050 FORMAT (' vacuum parameter=',f9.5,' change: ',l1,& - & ' second vacuum=',f9.5) - ENDIF - - RETURN - END SUBROUTINE w_enpara -! -!------------------------------------------------------------------ - SUBROUTINE r_enpara(& - & atoms,input,jsp,& - & enpara) -!------------------------------------------------------------------ - USE m_types - IMPLICIT NONE - - INTEGER, INTENT (IN) :: jsp - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_input),INTENT(IN) :: input - TYPE(t_enpara),INTENT(INOUT):: enpara - - INTEGER n,l,lo,skip_t,io_err - enpara%lchange(:,:,jsp) =.false. - enpara%el0(:,:,jsp) =0.0 - enpara%ello0(:,:,jsp) =0.0 - -!--> first line contains mixing parameter! - - enpara%enmix(jsp) = 0.0 - READ (40,FMT ='(48x,f10.6)',iostat=io_err) enpara%enmix(jsp) - IF (io_err /= 0) THEN - !use defaults - enpara%lchange(:,:,jsp)=.false. - enpara%llochg(:,:,jsp)=.false. - enpara%evac0(:,jsp) = eVac0Default_const - enpara%skiplo(:,jsp) = 0 - enpara%enmix(jsp) = 0.0 - enpara%lchg_v(:,jsp)=.false. - enpara%el0(:,:,jsp)=-99999. - CALL default_enpara(jsp,atoms,enpara) - WRITE (6,FMT=8001) jsp - WRITE (6,FMT = 8000) - DO n = 1,atoms%ntype - WRITE (6,FMT=8140) n,(enpara%el0(l,n,jsp),l=0,3),& - & (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp) - IF (atoms%nlo(n)>=1) THEN - WRITE (6,FMT = 8139) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n)) - WRITE (6,FMT = 8138) (enpara%llochg(lo,n,jsp),lo = 1,atoms%nlo(n)) - ENDIF - ENDDO - RETURN - ENDIF - READ (40,*) ! skip next line - IF (enpara%enmix(jsp).EQ.0.0) enpara%enmix(jsp) = 1.0 - WRITE (6,FMT=8001) jsp - WRITE (6,FMT=8000) - skip_t = 0 - DO n = 1,atoms%ntype - READ (40,FMT=8040,END=200) (enpara%el0(l,n,jsp),l=0,3),& - & (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp) - WRITE (6,FMT=8140) n,(enpara%el0(l,n,jsp),l=0,3),& - & (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp) -! -!---> energy parameters for the local orbitals -! - IF (atoms%nlo(n).GE.1) THEN - skip_t = skip_t + enpara%skiplo(n,jsp) * atoms%neq(n) - READ (40,FMT=8039,END=200) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n)) - READ (40,FMT=8038,END=200) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n)) - WRITE (6,FMT=8139) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n)) - WRITE (6,FMT=8138) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n)) - ELSEIF (enpara%skiplo(n,jsp).GT.0) THEN - WRITE (6,*) "for atom",n," no LO's were specified" - WRITE (6,*) 'but skiplo was set to',enpara%skiplo - CALL juDFT_error("No LO's but skiplo",calledby ="enpara",& - & hint="If no LO's are set skiplo must be 0 in enpara") - END IF -! -!---> set the energy parameters with l>3 to the value of l=3 -! - DO l = 4,atoms%lmax(n) - enpara%el0(l,n,jsp) = enpara%el0(3,n,jsp) - ENDDO - ENDDO ! atoms%ntype - - IF (input%film) THEN - enpara%lchg_v = .true. - READ (40,FMT=8050,END=200) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp) - WRITE (6,FMT=8150) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp) - ENDIF - IF (atoms%nlod.GE.1) THEN - WRITE (6,FMT=8090) jsp,skip_t - WRITE (6,FMT=8091) - END IF - -! input formats - - 8038 FORMAT (14x,60(l1,8x)) - 8039 FORMAT (8x,60f9.5) - 8040 FORMAT (8x,4f9.5,9x,4l1,9x,i3) - 8050 FORMAT (19x,f9.5,9x,l1,15x,f9.5) - -! output formats - - 8138 FORMAT (' --> change ',60(l1,8x)) - 8139 FORMAT (' --> lo ',60f9.5) - 8140 FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3) - 8150 FORMAT (' vacuum parameter=',f9.5,' change: ',l1,& - & ' second vacuum=',f9.5) - 8001 FORMAT ('READING enpara for spin: ',i1) - 8000 FORMAT (/,' energy parameters:',/,t10,'s',t20,& - & 'p',t30,'d',t37,'higher l - - -') - 8090 FORMAT ('Spin: ',i1,' -- ',i3,'eigenvalues') - 8091 FORMAT ('will be skipped for energyparameter computation') - - RETURN - - 200 WRITE (6,*) 'the end of the file enpara has been reached while' - WRITE (6,*) 'reading the energy-parameters.' - WRITE (6,*) 'possible reason: energy parameters have not been' - WRITE (6,*) 'specified for all atom types.' - WRITE (6,FMT='(a,i4)')& - & 'the actual number of atom-types is: ntype=',atoms%ntype - CALL juDFT_error& - & ("unexpected end of file enpara reached while reading"& - & ,calledby ="enpara") - - END SUBROUTINE r_enpara -! -!------------------------------------------------------------------ - SUBROUTINE mix_enpara(& - & jsp,atoms,vacuum,& - & obsolete,input,enpara,& - & vr,vz,pvac,svac,& - & ener,sqal,enerlo,sqlo) -!------------------------------------------------------------------ - USE m_types - IMPLICIT NONE - INTEGER,INTENT(IN) :: jsp - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_vacuum),INTENT(IN) :: vacuum - TYPE(t_obsolete),INTENT(IN) :: obsolete - TYPE(t_input),INTENT(IN) :: input - TYPE(t_enpara),INTENT(INOUT) :: enpara - - REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype) - REAL, INTENT(IN) :: ener(0:3,atoms%ntype),sqal(0:3,atoms%ntype) - REAL, INTENT(IN) :: enerlo(atoms%nlod,atoms%ntype),sqlo(atoms%nlod,atoms%ntype) - REAL, INTENT(IN) :: pvac(2),svac(2),vz(vacuum%nmzd,2) - - INTEGER ityp,j,l,lo - REAl vbar,maxdist - INTEGER same(atoms%nlod) - LOGICAL l_int_enpara - - l_int_enpara=all(enpara%el0==int(enpara%el0)) !test only enpara of lapw, no lo - - maxdist=0.0 - DO ityp = 1,atoms%ntype -! look for LO's energy parameters equal to the LAPW (and previous LO) ones - same = 0 - DO lo = 1,atoms%nlo(ityp) - IF(enpara%el0(atoms%llo(lo,ityp),ityp,jsp).eq.enpara%ello0(lo,ityp,jsp)) same(lo)=-1 - DO l = 1,lo-1 - IF(atoms%llo(l,ityp).ne.atoms%llo(lo,ityp)) cycle - IF(enpara%ello0(l,ityp,jsp).eq.enpara%ello0(lo,ityp,jsp).and.same(lo).eq.0)& - & same(lo)=l - ENDDO - ENDDO -! -!---> change energy parameters -! - IF ( obsolete%lepr.EQ.1) THEN - j = atoms%jri(ityp) - (log(4.0)/atoms%dx(ityp)+1.51) - vbar = vr(j,ityp)/( atoms%rmt(ityp)*exp(atoms%dx(ityp)*(j-atoms%jri(ityp))) ) - ELSE - vbar = 0.0 - END IF - - DO l = 0,3 - IF ( enpara%lchange(l,ityp,jsp) ) THEN - write(6,*) 'Type:',ityp,' l:',l - write(6,FMT=777) enpara%el0(l,ityp,jsp),& - & (ener(l,ityp)/sqal(l,ityp) - vbar),& - & abs(enpara%el0(l,ityp,jsp)-(ener(l,ityp)/sqal(l,ityp) - vbar)) - maxdist=max(maxdist,& - & abs(enpara%el0(l,ityp,jsp)-(ener(l,ityp)/sqal(l,ityp) - vbar))) - enpara%el0(l,ityp,jsp) =(1.0-enpara%enmix(jsp))*enpara%el0(l,ityp,jsp) + & - & enpara%enmix(jsp)*(ener(l,ityp)/sqal(l,ityp) - vbar) - ENDIF - ENDDO - DO l = 4, atoms%lmaxd - IF ( enpara%lchange(3,ityp,jsp) ) THEN - enpara%el0(l,ityp,jsp) = enpara%el0(3,ityp,jsp) - ENDIF - ENDDO -! -!---> determine and change local orbital energy parameters -! - DO lo = 1,atoms%nlo(ityp) - IF (atoms%l_dulo(lo,ityp)) THEN - enpara%ello0(lo,ityp,jsp) =enpara%el0(atoms%llo(lo,ityp),ityp,jsp) - ELSEIF (enpara%llochg(lo,ityp,jsp) ) THEN - IF(same(lo).eq.-1) THEN - enpara%ello0(lo,ityp,jsp) = enpara%el0(atoms%llo(lo,ityp),ityp,jsp) - cycle - ELSE IF(same(lo).gt.0) THEN - enpara%ello0(lo,ityp,jsp) = enpara%ello0(same(lo),ityp,jsp) - cycle - ENDIF - write(6,*) 'Type:',ityp,' lo:',lo - write(6,FMT=777) enpara%ello0(lo,ityp,jsp),& - & (enerlo(lo,ityp)/sqlo(lo,ityp) - vbar),& - & abs(enpara%ello0(lo,ityp,jsp)-(enerlo(lo,ityp)/sqlo(lo,ityp)-vbar)) - maxdist=max(maxdist,& - & abs(enpara%ello0(lo,ityp,jsp)-(enerlo(lo,ityp)/sqlo(lo,ityp)-vbar))) - enpara%ello0(lo,ityp,jsp) =(1.0-enpara%enmix(jsp))*enpara%ello0(lo,ityp,jsp)+& - & enpara%enmix(jsp)*(enerlo(lo,ityp)/sqlo(lo,ityp) - vbar) - END IF - END DO - - END DO - - IF (input%film) THEN - - IF (obsolete%lepr.eq.1) THEN - vbar = vz(1,1) - ELSE - vbar = 0.0 - ENDIF - - IF (enpara%lchg_v(1,jsp) ) THEN - write(6,*) 'Vacuum:' - write(6,FMT=777) enpara%evac0(1,jsp),(pvac(1)/svac(1) - vbar),& - & abs(enpara%evac0(1,jsp)-(pvac(1)/svac(1) - vbar)) - maxdist=max(maxdist,abs(enpara%evac0(1,jsp)-(pvac(1)/svac(1) - vbar))) - enpara%evac0(1,jsp) =(1.0-enpara%enmix(jsp))*enpara%evac0(1,jsp)+& - & enpara%enmix(jsp)*(pvac(1)/svac(1) - vbar) - IF (vacuum%nvac.EQ.2) THEN - IF (obsolete%lepr.eq.1) vbar = vz(1,vacuum%nvac) - write(6,FMT=777) enpara%evac0(2,jsp),(pvac(2)/svac(2) - vbar),& - & abs(enpara%evac0(2,jsp)-(pvac(2)/svac(2) - vbar)) - maxdist=max(maxdist,abs(enpara%evac0(2,jsp)-(pvac(2)/svac(2) - vbar))) - enpara%evac0(2,jsp) = (1.0-enpara%enmix(jsp))*enpara%evac0(2,jsp)+& - & enpara%enmix(jsp)*(pvac(2)/svac(2) - vbar) - ELSE - enpara%evac0(2,jsp) = enpara%evac0(1,jsp) - ENDIF - ENDIF - - ENDIF - WRITE(6,'(a36,f12.6)') 'Max. mismatch of energy parameters:',& - & maxdist - if (maxdist>1.0.and..not.l_int_enpara) call juDFT_warn& - & ("Energy parameter mismatch too large",hint& - & ="If any energy parameters calculated from the output "//& - & "differ from the input by more than 1Htr, chances are "//& - & "high that your initial setup was broken.") - RETURN - 777 FORMAT('Old:',f8.5,' new:',f8.5,' diff:',f8.5) - - END SUBROUTINE mix_enpara - - subroutine default_enpara(jsp,atoms,enpara) - !Create default values for enpara given as simple integers - use m_types - implicit none - INTEGER ,INTENT(IN) :: jsp - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_enpara),INTENT(INOUT):: enpara - - INTEGER :: n,i - - - DO n = 1,atoms%ntype - IF (all(enpara%el0(:,n,jsp)>-9999.)) cycle !enpara was set already - IF ( atoms%nz(n) < 3 ) THEN - enpara%el0(:,n,jsp) = real( (/1,2,3,4/) ) - ELSEIF ( atoms%nz(n) < 11 ) THEN - enpara%el0(:,n,jsp) = real( (/2,2,3,4/) ) - ELSEIF ( atoms%nz(n) < 19 ) THEN - enpara%el0(:,n,jsp) = real( (/3,3,3,4/) ) - ELSEIF ( atoms%nz(n) < 31 ) THEN - enpara%el0(:,n,jsp) = real( (/4,4,3,4/) ) - ELSEIF ( atoms%nz(n) < 37 ) THEN - enpara%el0(:,n,jsp) = real( (/4,4,4,4/) ) - ELSEIF ( atoms%nz(n) < 49 ) THEN - enpara%el0(:,n,jsp) = real( (/5,5,4,4/) ) - ELSEIF ( atoms%nz(n) < 55 ) THEN - enpara%el0(:,n,jsp) = real( (/5,5,5,4/) ) - ELSEIF ( atoms%nz(n) < 72 ) THEN - enpara%el0(:,n,jsp) = real( (/6,6,5,4/) ) - ELSEIF ( atoms%nz(n) < 81 ) THEN - enpara%el0(:,n,jsp) = real( (/6,6,5,5/) ) - ELSEIF ( atoms%nz(n) < 87 ) THEN - enpara%el0(:,n,jsp) = real( (/6,6,6,5/) ) - ELSE - enpara%el0(:,n,jsp) = real( (/7,7,6,5/) ) - ENDIF - DO i = 1, atoms%nlo(n) - enpara%ello0(i,n,jsp) = enpara%el0(atoms%llo(i,n),n,jsp) - 1 - ENDDO - ENDDO - END subroutine - END MODULE m_enpara diff --git a/global/nvtx.F90 b/global/nvtx.F90 deleted file mode 100644 index 3fc476043c7d1bb99cb8f8808c49757584e53839..0000000000000000000000000000000000000000 --- a/global/nvtx.F90 +++ /dev/null @@ -1,64 +0,0 @@ -module nvtx -! See https://devblogs.nvidia.com/parallelforall/customize-cuda-fortran-profiling-nvtx/ -use iso_c_binding -implicit none - -integer,private :: col(7) = [ Z'0000ff00', Z'000000ff', Z'00ffff00', Z'00ff00ff', Z'0000ffff', Z'00ff0000', Z'00ffffff'] -character(len=256),private :: tempName - -type, bind(C):: nvtxEventAttributes - integer(C_INT16_T):: version=1 - integer(C_INT16_T):: size=48 ! - integer(C_INT):: category=0 - integer(C_INT):: colorType=1 ! NVTX_COLOR_ARGB = 1 - integer(C_INT):: color - integer(C_INT):: payloadType=0 ! NVTX_PAYLOAD_UNKNOWN = 0 - integer(C_INT):: reserved0 - integer(C_INT64_T):: payload ! union uint,int,double - integer(C_INT):: messageType=1 ! NVTX_MESSAGE_TYPE_ASCII = 1 - type(C_PTR):: message ! ascii char -end type - -interface nvtxRangePush - ! push range with custom label and standard color - subroutine nvtxRangePushA(name) bind(C, name='nvtxRangePushA') - use iso_c_binding - character(kind=C_CHAR,len=*) :: name - end subroutine - - ! push range with custom label and custom color - subroutine nvtxRangePushEx(event) bind(C, name='nvtxRangePushEx') - use iso_c_binding - import:: nvtxEventAttributes - type(nvtxEventAttributes):: event - end subroutine -end interface - -interface nvtxRangePop - subroutine nvtxRangePop() bind(C, name='nvtxRangePop') - end subroutine -end interface - -contains - -subroutine nvtxStartRange(name,id) - character(kind=c_char,len=*) :: name - integer, optional:: id - type(nvtxEventAttributes):: event - - tempName=trim(name)//c_null_char - - if ( .not. present(id)) then - call nvtxRangePush(tempName) - else - event%color=col(mod(id,7)+1) - event%message=c_loc(tempName) - call nvtxRangePushEx(event) - end if -end subroutine - -subroutine nvtxEndRange - call nvtxRangePop -end subroutine - -end module nvtx diff --git a/io/eig66_parallel.F90 b/io/eig66_parallel.F90 deleted file mode 100644 index 073632ca6f32019f260fc0f3ba1b2acbf755a2a2..0000000000000000000000000000000000000000 --- a/io/eig66_parallel.F90 +++ /dev/null @@ -1,200 +0,0 @@ -module m_eig66_parallel -#include "juDFT_env.h" -! Do the IO of the eig-file into memory, different from simple approach by -! using k-point parallelism -! The eig-file is split into four arrays: -! eig_int contains the basis-set information/integers (nv,nmat,ne,k1,k2,k3,kveclo) -! eig_real contains the basis-set information/real (el,evac,ello,bkpt,wtkpt) -! eig_eig contains the eigenvalues -! eig_vec contains the eigenvectors -! The record number is given by nrec=nk+(jspin-1)*nkpts - implicit none - INTEGER,ALLOCATABLE :: eig_int(:,:) - REAL,ALLOCATABLE :: eig_real(:,:) - REAL,ALLOCATABLE :: eig_eig(:,:) - REAL,ALLOCATABLE :: eig_vecr(:,:) - COMPLEX,ALLOCATABLE :: eig_vecc(:,:) - integer :: size_k,isize - - INTEGER:: ntypes,lmax,nlo,jspins,nkpts - - contains - subroutine open_eig(nmat,neig,nkpts_in,jspins_in,lmax_in,nlo_in,ntype_in,l_create,nlotot,l_noco) - INTEGER, INTENT(IN) :: nmat,neig,nkpts_in,jspins_in,nlo_in,ntype_in,lmax_in,nlotot - LOGICAL, INTENT(IN) :: l_noco,l_create - - integer:: length,e - - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,isize,e) - - if (mod(nkpts_in/isize).ne.0) call judft_error("Could not distribute k-points in eig66_parallel") - - if (allocated(eig_int)) then - if (.not.l_create) return - call close_eig(.true.) - endif - - - ntypes=ntype_in - lmax=lmax_in - nlo=nlo_in - jspins=jspins_in - nkpts=nkpts_in/isize - - !eig_int - length=3 !nv+nmat+ne - length=length+nmat*3 !k1,k2,k3 - length=length+nlotot !kveclo - size_k=nmat - allocate(eig_int(length,jspins*nkpts)) - - !eig_real - length=3+1+2 !bk,wk,evac - length=length+lmax*ntypes !el - length=length+nlo*ntypes !ello - ALLOCATE(eig_real(length,jspins*nkpts)) - !eig_eig - length=jspins - if (l_noco) length=1 - allocate(eig_eig(neig,length*nkpts)) - !eig_vec -#ifdef CPP_INVERSION - allocate(eig_vecr(nmat*neig,length*nkpts)) -#ifdef CPP_SOC - call judft_error("SOC+INVERSION can not be used with eigenvalues stored in memory") -#endif -#else - allocate(eig_vecc(nmat*neig,length*nkpts)) -#endif - end subroutine open_eig - - subroutine close_eig(delete) - logical,intent(in),optional::delete - if (present(delete)) THEN - if (delete) THEN - if (allocated(eig_int)) deallocate(eig_int) - if (allocated(eig_real)) deallocate(eig_real) - if (allocated(eig_eig)) deallocate(eig_eig) - if (allocated(eig_vecr)) deallocate(eig_vecr) - if (allocated(eig_vecc)) deallocate(eig_vecc) - endif - endif - - end subroutine close_eig - - subroutine read_eig(nk,jspin,nv,nmat,k1,k2,k3,bk,wk,neig,eig,el,& - ello,evac,kveclo,n_start,n_end,zr,zc) - implicit none - INTEGER, INTENT(IN) :: nk,jspin - INTEGER, INTENT(OUT),OPTIONAL :: nv,nmat - INTEGER, INTENT(OUT),OPTIONAL :: neig - REAL, INTENT(OUT),OPTIONAL :: eig(:) - INTEGER, INTENT(OUT),OPTIONAL :: k1(:),k2(:),k3(:),kveclo(:) - REAL, INTENT(OUT),OPTIONAL :: evac(:),ello(:,:),el(:,:) - REAL, INTENT(OUT),OPTIONAL :: bk(:),wk - INTEGER, INTENT(IN),OPTIONAL :: n_start,n_end - REAL, INTENT(OUT),OPTIONAL :: zr(:,:) - COMPLEX, INTENT(OUT),OPTIONAL :: zc(:,:) - - INTEGER::nrec - - nrec=nk/isize+(jspin-1)*nkpts/isize+1 - ! data from eig_int - if (present(nv)) nv=eig_int(1,nrec) - if (present(nmat)) nmat=eig_int(2,nrec) - if (present(neig)) then - neig=eig_int(3,nrec) - endif - if (present(k1)) then - if (.not.present(k2).or..not.present(k3)) call juDFT_error("BUG: always read k1,k2,k3") - k1=eig_int(3+1:3+size_k,nrec) - k2=eig_int(3+size_k+1:3+2*size_k,nrec) - k3=eig_int(3+2*size_k+1:3+3*size_k,nrec) - endif - if (present(kveclo)) kveclo=eig_int(4+3*size_k:3+3*size_k+size(kveclo),nrec) - - !data from eig_real - if (present(bk)) bk=eig_real(1:3,nrec) - if (present(wk)) wk=eig_real(4,nrec) - if (present(evac)) evac=eig_real(5:6,nrec) - if (present(el)) el=reshape(eig_real(7:7+size(el)-1,nrec),shape(el)) - if (present(ello)) ello=reshape(eig_real(size(eig_real,1)-size(ello)+1:,nrec),shape(ello)) - - !data from eig_eig - if (present(eig)) THEN - eig=0.0 - eig=eig_eig(:size(eig),nrec) - ENDIF - !data from eig_vec - if (present(zr)) then - write(*,*) "R-Z:",nrec,shape(zr) - - if (.not.allocated(eig_vecr)) call juDFT_error("BUG: can not read real vectors from memory") - zr=reshape(eig_vecr(:size(zr),nrec),shape(zr)) - - endif - if (present(zc)) then - write(*,*) "R-ZC:",nrec,shape(zc) - - if (.not.allocated(eig_vecc)) call juDFT_error("BUG: can not read complex vectors from memory") - zc=reshape(eig_vecc(:size(zc),nrec),shape(zc)) - endif - end subroutine read_eig - - - subroutine write_eig(nk,jspin,neig,nv,nmat,k1,k2,k3,bk,wk, & - eig,el,ello,evac, & - nlotot,kveclo,n_size,n_rank,zr,zc) - INTEGER, INTENT(IN) :: nk,jspin - INTEGER, INTENT(IN),OPTIONAL :: n_size,n_rank - REAL, INTENT(IN),OPTIONAL :: wk - INTEGER, INTENT(IN),OPTIONAL :: neig,nv,nmat,nlotot - INTEGER, INTENT(IN),OPTIONAL :: k1(:),k2(:),k3(:),kveclo(:) - REAL, INTENT(IN),OPTIONAL :: bk(3),eig(:),el(:,:) - REAL, INTENT(IN),OPTIONAL :: evac(:),ello(:,:) - REAL, INTENT(IN),OPTIONAL :: zr(:,:) - COMPLEX,INTENT(IN),OPTIONAL :: zc(:,:) - INTEGER::nrec - - nrec=(nk+(jspin-1)*nkpts)/isize+1 - ! data from eig_int - if (present(nv)) eig_int(1,nrec)=nv - if (present(nmat)) eig_int(2,nrec)=nmat - if (present(neig)) THEN - eig_int(3,nrec)=neig - endif - if (present(k1)) then - if (.not.present(k2).or..not.present(k3)) call juDFT_error("BUG: always write k1,k2,k3") - eig_int(3+1:3+size_k,nrec)=k1 - eig_int(3+size_k+1:3+2*size_k,nrec)=k2 - eig_int(3+2*size_k+1:3+3*size_k,nrec)=k3 - endif - if (present(kveclo)) eig_int(4+3*size_k:3+3*size_k+size(kveclo),nrec)=kveclo - - !data from eig_real - if (present(bk)) eig_real(1:3,nrec)=bk - if (present(wk)) eig_real(4,nrec)=wk - if (present(evac)) eig_real(5:6,nrec)=evac - if (present(el)) eig_real(7:7+size(el)-1,nrec)=reshape(el,(/size(el)/)) - if (present(ello)) eig_real(size(eig_real,1)-size(ello)+1:,nrec)=reshape(ello,(/size(ello)/)) - !data from eig_eig - if (present(eig)) eig_eig(:size(eig),nrec)=eig - !data from eig_vec - if (present(zr)) then - write(*,*) "W-ZR:",nrec,shape(zr) - - if (.not.allocated(eig_vecr)) call juDFT_error("BUG: can not write real vectors to memory") - eig_vecr(:size(zr),nrec)=reshape(zr,(/size(zr)/)) - endif - if (present(zc)) then - write(*,*) "W-ZC:",nrec,shape(zc) - - if (.not.allocated(eig_vecc)) call juDFT_error("BUG: can not write complex vectors to memory") - eig_vecc(:size(zc),nrec)=reshape(zc,(/size(zc)/)) - endif - - - end subroutine write_eig - - -end module m_eig66_mem diff --git a/optional/f2u.f90 b/optional/f2u.f90 deleted file mode 100644 index c05ba8c6dd1ab1ddc384e7c99aa577a52378cedb..0000000000000000000000000000000000000000 --- a/optional/f2u.f90 +++ /dev/null @@ -1,165 +0,0 @@ - MODULE m_f2u - use m_juDFT -! **************************************************** -! write unformatted density or potential onto unit '99' -! e. wimmer march 1985 -! **************************************************** - CONTAINS - SUBROUTINE f2u(& - & stars,input,atoms,sphhar,vacuum,& - & cell,sym,l_noco) -! - USE m_wrtdop - USE m_types - IMPLICIT NONE -! .. -! .. Scalar Arguments .. - TYPE(t_stars),INTENT(IN) :: stars - TYPE(t_input),INTENT(IN) :: input - TYPE(t_vacuum),INTENT(IN) :: vacuum - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_sphhar),INTENT(IN) :: sphhar - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_cell),INTENT(IN) :: cell - LOGICAL,INTENT(IN) :: l_noco - -! .. -! -! .. Local Scalars .. - INTEGER it,nu,i,jsp,jspdum,nn,n,ndum,jrin,j,ntypsyn,na - INTEGER nlhn,lh,lhdum,nq3n,k,ivac,ivdum,nmzn,nq2n,nmzxyn - REAL z1n,delzn,rmtn,dxn - LOGICAL n_exist - CHARACTER(len=2) namaux - CHARACTER(len=8) dop,iop -! .. -! .. Local Arrays .. - COMPLEX, ALLOCATABLE :: n_mmp(:,:,:) - COMPLEX, ALLOCATABLE :: fpw(:,:),fzxy(:,:,:,:) - COMPLEX, ALLOCATABLE :: cdom(:),cdomvz(:,:),cdomvxy(:,:,:) - REAL, ALLOCATABLE :: fz(:,:,:),fr(:,:,:,:) - REAL, ALLOCATABLE :: fpwr(:,:),fzxyr(:,:,:,:) - CHARACTER(len=8) name(10),space(10) -! .. -! .. Data statements .. - DATA space/10*' '/ -! .. -! - !atoms%jmtd = maxval(atoms%jri(:)) - !sphhar%nlhd = maxval(sphhar%nlh(:)) - - IF (l_noco) THEN - OPEN (98,file='rhomat_inp',form='unformatted') - OPEN (99,file='rhomat_form',form='formatted') - ELSE - OPEN (98,file='cdn_unf',form='unformatted') - OPEN (99,file='cdn_form',form='formatted') - ENDIF - nu=99 -! - ALLOCATE( fpw(stars%ng3,input%jspins),fzxy(vacuum%nmzxy,stars%ng2-1,2,input%jspins) ) - ALLOCATE( cdom(stars%ng3),cdomvz(vacuum%nmz,2),cdomvxy(vacuum%nmzxy,stars%ng2-1,2) ) - ALLOCATE( fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins),fz(vacuum%nmz,2,input%jspins) ) - ALLOCATE( fpwr(stars%ng3,input%jspins),fzxyr(vacuum%nmzxy,stars%ng2-1,2,input%jspins) ) - - DO 10 i = 1,10 - name(i) = space(i) - 10 CONTINUE - READ (nu,FMT=8030) name - WRITE (6,FMT=8000) name - 8000 FORMAT (' loddop title:',10a8) - READ (nu,FMT=8050) iop,dop,it - DO 60 jsp = 1,input%jspins - READ (nu,FMT=8060) jspdum - READ (nu,FMT=8070) nn - IF (nn/=atoms%ntype) CALL juDFT_error("loddop",calledby ="f2u") - na = 1 - DO 30 n = 1,nn - READ (nu,FMT=8010) namaux,ndum,jrin,rmtn,dxn - 8010 FORMAT (a2,4x,i3,6x,i3,6x,f10.6,5x,f10.6) - READ (nu,FMT=8020) ntypsyn,nlhn - IF (ntypsyn.NE.atoms%ntypsy(na)) THEN - WRITE (6,*) 'ntypsyn = ',ntypsyn, ' =/=' - WRITE (6,*) 'ntypsy(na) = ',atoms%ntypsy(na) - CALL juDFT_error("ntypsyn /= ntypsy(na)",calledby ="f2u") - ENDIF - 8020 FORMAT (7x,i2,2x,4x,i3) - DO 20 lh = 0,nlhn - READ (nu,FMT=8060) lhdum - READ (nu,FMT=8100) (fr(i,lh,n,jsp),i=1,jrin) - 20 CONTINUE - na = na + atoms%neq(n) - 30 CONTINUE - READ (nu,FMT=8110) nq3n - IF (sym%invs) THEN - READ (nu,FMT=8100) (fpwr(k,jsp),k=1,nq3n) - DO k=1,nq3n - fpw(k,jsp) = cmplx(fpwr(k,jsp),0.0) - ENDDO - ELSE - READ (nu,FMT=8100) (fpw(k,jsp),k=1,nq3n) - END IF - IF (input%film) THEN - DO 50 ivac = 1,vacuum%nvac - READ (nu,FMT=8060) ivdum - READ (nu,FMT=8120) nmzn,z1n,delzn - READ (nu,FMT=8040) (fz(i,ivac,jsp),i=1,nmzn) - READ (nu,FMT=8130) nq2n,nmzxyn - DO 40 k = 2,nq2n - IF (sym%invs2) THEN - READ (nu,FMT=8100) (fzxyr(j,k-1,ivac,jsp),j=1,& - & nmzxyn) - DO j=1,nmzxyn - fzxy(j,k-1,ivac,jsp) = & - & cmplx(fzxyr(j,k-1,ivac,jsp),0.0) - ENDDO - ELSE - READ (nu,FMT=8100) (fzxy(j,k-1,ivac,jsp),j=1,& - & nmzxyn) - END IF - 40 CONTINUE - 50 CONTINUE - END IF - 60 CONTINUE - - IF (l_noco) THEN - READ (99,8100) (cdom(k),k=1,stars%ng3) - IF (input%film) THEN - READ (99,8100) ((cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac) - READ (99,8100) (((cdomvxy(j,k-1,ivac),j=1,vacuum%nmzxy),k=2,stars%ng2),& - & ivac=1,vacuum%nvac) - ENDIF - ENDIF - - CLOSE (99) - - 8030 FORMAT (10a8) - 8040 FORMAT (4e20.13) - 8050 FORMAT (a8,1x,a8,11x,i3) - 8060 FORMAT (a9) - 8070 FORMAT (6x,i3) - 8080 FORMAT (15x,i3,6x,f10.6,5x,f10.6) - 8090 FORMAT (15x,i2) - 8100 FORMAT (4e20.13) - 8110 FORMAT (4x,i6) - 8120 FORMAT (4x,i3,5x,f20.13,7x,f8.4) - 8130 FORMAT (4x,i5,8x,i5) -! - CALL wrtdop(& - & stars,vacuum,atoms,sphhar,& - & input,sym,98,& - & it,fr,fpw,fz,fzxy) - - IF (l_noco) THEN - WRITE (98) (cdom(k),k=1,stars%ng3) - IF (input%film) THEN - WRITE (98) ((cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac) - WRITE (98) (((cdomvxy(j,k-1,ivac),j=1,vacuum%nmzxy),k=2,stars%ng2),& - & ivac=1,vacuum%nvac) - ENDIF - ENDIF - CLOSE (98) - DEALLOCATE( fpw,fzxy,cdom,cdomvz,cdomvxy,fr,fz,fpwr,fzxyr ) - - END SUBROUTINE f2u - END MODULE m_f2u diff --git a/optional/u2f.f90 b/optional/u2f.f90 deleted file mode 100644 index 4f94529b7f497ab778801dcc89018bff0802103f..0000000000000000000000000000000000000000 --- a/optional/u2f.f90 +++ /dev/null @@ -1,141 +0,0 @@ - MODULE m_u2f -! **************************************************** -! write formatted density or potential onto unit '99' -! e. wimmer march 1985 -! **************************************************** - CONTAINS - SUBROUTINE u2f(& - & stars,input,atoms,sphhar,vacuum,& - & cell,sym,l_noco) -! - USE m_loddop - USE m_constants,ONLY:namat_const - USE m_types - IMPLICIT NONE -! .. -! .. Scalar Arguments .. - TYPE(t_stars),INTENT(IN) :: stars - TYPE(t_input),INTENT(IN) :: input - TYPE(t_vacuum),INTENT(IN) :: vacuum - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_sphhar),INTENT(IN) :: sphhar - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_cell),INTENT(IN) :: cell - LOGICAL,INTENT(IN) :: l_noco -! .. -! .. -! .. Local Scalars .. - INTEGER i,ivac,j,jsp,k,lh,n,izn,na,urec - LOGICAL n_exist -! .. -! .. Local Arrays .. - INTEGER it - CHARACTER(len=8) dop,iop,name(10) - COMPLEX, ALLOCATABLE :: fpw(:,:),fzxy(:,:,:,:) - COMPLEX, ALLOCATABLE :: cdom(:),cdomvz(:,:),cdomvxy(:,:,:) - REAL, ALLOCATABLE :: fz(:,:,:),fr(:,:,:,:) -! .. -! .. - !atoms%jmtd = maxval(atoms%jri(:)) - !sphhar%nlhd = maxval(sphhar%nlh(:)) - - IF (l_noco) THEN - OPEN (98,file='rhomat_inp',form='unformatted') - OPEN (99,file='rhomat_form',form='formatted') - ELSE - OPEN (98,file='f_unf',form='unformatted') - OPEN (99,file='f_form',form='formatted') - ENDIF -! - ALLOCATE( fpw(stars%ng3,input%jspins),fzxy(vacuum%nmzxy,stars%ng2-1,2,input%jspins) ) - ALLOCATE( cdom(stars%ng3),cdomvz(vacuum%nmz,2),cdomvxy(vacuum%nmzxy,stars%ng2-1,2) ) - ALLOCATE( fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins),fz(vacuum%nmz,2,input%jspins) ) - - CALL loddop(& - & stars,vacuum,atoms,sphhar,& - & input,sym,& - & 98,& - & it,fr,fpw,fz,fzxy) - write(*,*) 'after loddop:', l_noco - IF (l_noco) THEN - cdom(:) = (0.0,0.0) - cdomvz(:,:) = (0.0,0.0) ; cdomvxy(:,:,:) = (0.0,0.0) - READ (98,END=101,ERR=101) (cdom(k),k=1,stars%ng3) - IF (input%film) THEN - READ (98) ((cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac) - READ (98) (((cdomvxy(j,k-1,ivac),j=1,vacuum%nmzxy),k=2,stars%ng2),& - & ivac=1,vacuum%nvac) - ENDIF - ENDIF - - 101 CLOSE (98) - - WRITE (99,FMT=8000) name - 8000 FORMAT (10a8) - WRITE (6,FMT=8010) name - 8010 FORMAT (' loddop title:',10a8) - WRITE (99,FMT=8020) iop,dop,it - 8020 FORMAT (a8,1x,a8,' iteration=',i3) - DO 60 jsp = 1,input%jspins - WRITE (99,FMT=8030) jsp - 8030 FORMAT ('spin=',i1) - WRITE (99,FMT=8040) atoms%ntype - 8040 FORMAT ('ntype=',i3) - na = 1 - DO 30 n = 1,atoms%ntype - izn = atoms%zatom(n) + 0.01 - WRITE (99,FMT=8050) namat_const(izn),n,atoms%jri(n),atoms%rmt(n),atoms%dx(n) - 8050 FORMAT (a2,2x,'n=',i3,2x,'jri=',i3,2x,'rmt=',f10.6,2x,'dx=',& - & f10.6) - WRITE (99,FMT=8060) atoms%ntypsy(na),sphhar%nlh(atoms%ntypsy(na)) - 8060 FORMAT ('ntypsy=',i2,2x,'nlh=',i3) - DO 20 lh = 0,sphhar%nlh(atoms%ntypsy(na)) - WRITE (99,FMT=8070) lh - 8070 FORMAT ('lh=',i3) - WRITE (99,FMT=8080) (fr(i,lh,n,jsp),i=1,atoms%jri(n)) - 8080 FORMAT (4e20.13) - 20 CONTINUE - na = na + atoms%neq(n) - 30 CONTINUE - WRITE (99,FMT=8090) stars%ng3 - 8090 FORMAT ('nq3=',i6) - IF (sym%invs) THEN - WRITE (99,FMT=8080) (real(fpw(k,jsp)),k=1,stars%ng3) - ELSE - WRITE (99,FMT=8080) (fpw(k,jsp),k=1,stars%ng3) - END IF - IF (input%film) THEN - DO 50 ivac = 1,vacuum%nvac - WRITE (99,FMT=8100) ivac - 8100 FORMAT ('ivac=',i1) - WRITE (99,FMT=8110) vacuum%nmz,cell%z1,vacuum%delz - 8110 FORMAT ('nmz=',i3,2x,'z1=',f20.13,2x,'delz=',f8.4) - WRITE (99,FMT=8080) (fz(i,ivac,jsp),i=1,vacuum%nmz) - WRITE (99,FMT=8120) stars%ng2,vacuum%nmzxy - 8120 FORMAT ('nq2=',i5,2x,'nmzxy=',i5) - DO 40 k = 2,stars%ng2 - IF (sym%invs2) THEN - WRITE (99,FMT=8080) (real(fzxy(j,k-1,ivac,jsp)),& - & j=1,vacuum%nmzxy) - ELSE - WRITE (99,FMT=8080) (fzxy(j,k-1,ivac,jsp),j=1,& - & vacuum%nmzxy) - END IF - 40 CONTINUE - 50 CONTINUE - END IF - 60 CONTINUE - - IF (l_noco) THEN - WRITE (99,8080) (cdom(k),k=1,stars%ng3) - IF (input%film) THEN - WRITE (99,8080) ((cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac) - WRITE (99,8080) (((cdomvxy(j,k-1,ivac),j=1,vacuum%nmzxy),k=2,stars%ng2),& - & ivac=1,vacuum%nvac) - ENDIF - ENDIF - CLOSE (99) - DEALLOCATE( fpw,fzxy,cdom,cdomvz,cdomvxy,fr,fz) - - END SUBROUTINE u2f - END MODULE m_u2f diff --git a/vgen/rhodirgen.f90 b/vgen/rhodirgen.f90 deleted file mode 100644 index bb37edfa2bd2283c4cf0dce7371ac1f8af903a49..0000000000000000000000000000000000000000 --- a/vgen/rhodirgen.f90 +++ /dev/null @@ -1,330 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- - -MODULE m_rhodirgen - USE m_juDFT - !********************************************************************** - ! This subroutine calculates the spin-up and -down density, which - ! are needed to calculate the potential and writes them to the file - ! cdn. The local angle of the magnetization is kept in real space - ! and written to the file dirofmag. This is done in four steps. - ! - ! i) The components of the hermitian density matrix (rho_11, rho_22, - ! rho_21) are reloaded from the file rhomat_inp. - ! ii) The density matrix in fouriertransformed to real space. - ! iii) The spin-up and -down densities and the local angle of the - ! magnetization are calculated on the real space mesh. - ! iv) The spin-up and -down densities are Fouriertransformed, stored - ! in terms of stars and written to the file cdn. The local angle of - ! magnetization is kept on the real space mesh and written to the - ! file dirofmag. - ! - ! Philipp Kurz 99/11/01 - !********************************************************************** -CONTAINS - SUBROUTINE rhodirgen(DIMENSION,sym,stars,atoms,sphhar,vacuum,& - cell,input,noco,oneD,den) - - !******** ABBREVIATIONS *********************************************** - ! ifft3 : size of the 3d real space mesh - ! ifft2 : size of the 2d real space mesh - ! rpw : diagonal components of the density matrix (rho_11 , - ! rho_22) - ! later interstitial spin-up and -down density - ! all stored in terms of 3d-stars - ! ris : first components of the density matrix - ! later interstitial spin-up and -down density and - ! direction of magnetic field (theta and phi) - ! all stored on real space mesh - !********************************************************************** - - USE m_constants - USE m_fft2d - USE m_fft3d - USE m_types - IMPLICIT NONE - - TYPE(t_dimension),INTENT(IN) :: DIMENSION - TYPE(t_noco),INTENT(IN) :: noco - TYPE(t_oneD),INTENT(IN) :: oneD - TYPE(t_input),INTENT(IN) :: input - TYPE(t_vacuum),INTENT(IN) :: vacuum - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_stars),INTENT(IN) :: stars - TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_sphhar),INTENT(IN) :: sphhar - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_potden),INTENT(INOUT) :: den - - ! .. Local Scalars .. - INTEGER iden,jspin,ivac,ifft2,ifft3 - INTEGER imz,ityp,iri,ilh,imesh,iq2,iq3 - REAL rho_11,rho_22,rho_21r,rho_21i,rhotot,magmom,phi - REAL rho_up,rho_down,mx,my,mz,eps,vz_r,vz_i,rziw,theta - ! .. - ! .. Local Arrays .. - REAL, ALLOCATABLE :: rz(:,:,:) - REAL, ALLOCATABLE :: rvacxy(:,:,:,:),ris(:,:),fftwork(:) - ! .. - eps = 1.0e-20 - - - - ! - !---> initialize arrays for the density matrix - ! - - ifft3 = 27*stars%mx1*stars%mx2*stars%mx3 - IF (input%film) THEN - ifft2 = 9*stars%mx1*stars%mx2 - IF (oneD%odi%d1) ifft2 = 9*stars%mx3*oneD%odi%M - ELSE - ifft2=0 - END IF - - IF (ALLOCATED(den%phi_pw)) THEN - DEALLOCATE(den%phi_pw,den%phi_vacz,den%phi_vacxy) - DEALLOCATE(den%theta_pw,den%theta_vacz,den%theta_vacxy) - ENDIF - ALLOCATE(den%phi_pw(ifft3),den%theta_pw(ifft3)) - ALLOCATE(den%phi_vacz(vacuum%nmzd,2),den%theta_vacz(vacuum%nmzd,2)) - ALLOCATE(den%phi_vacxy(ifft2,vacuum%nmzxyd,2),den%theta_vacxy(ifft2,vacuum%nmzxyd,2)) - - - ALLOCATE (ris(ifft3,4),fftwork(ifft3)) - !---> fouriertransform the diagonal part of the density matrix - !---> in the interstitial, den%pw, to real space (ris) - DO iden = 1,2 - CALL fft3d(ris(:,iden),fftwork,den%pw(:,iden),stars,+1) - ENDDO - !---> fouriertransform the off-diagonal part of the density matrix - CALL fft3d(ris(:,3),ris(:,4),den%pw(:,3),stars,+1) - - !test - ! DO iden=1,4 - ! write(*,*)'iden=',iden - ! write(*,8500)(ris(imesh,iden),imesh=0,ifft3-1) - ! enddo - !test - !---> calculate the charge and magnetization density on the - !---> real space mesh - DO imesh = 1,ifft3 - rho_11 = ris(imesh,1) - rho_22 = ris(imesh,2) - rho_21r = ris(imesh,3) - rho_21i = ris(imesh,4) - mx = 2*rho_21r - my = -2*rho_21i - mz = (rho_11-rho_22) - magmom = SQRT(mx**2 + my**2 + mz**2) - rhotot = rho_11 + rho_22 - rho_up = (rhotot + magmom)/2 - rho_down= (rhotot - magmom)/2 - - IF (ABS(mz) .LE. eps) THEN - theta = pi_const/2 - ELSEIF (mz .GE. 0.0) THEN - theta = ATAN(SQRT(mx**2 + my**2)/mz) - ELSE - theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const - ENDIF - - IF (ABS(mx) .LE. eps) THEN - IF (ABS(my) .LE. eps) THEN - phi = 0.0 - ELSEIF (my .GE. 0.0) THEN - phi = pi_const/2 - ELSE - phi = -pi_const/2 - ENDIF - ELSEIF (mx .GE. 0.0) THEN - phi = ATAN(my/mx) - ELSE - IF (my .GE. 0.0) THEN - phi = ATAN(my/mx) + pi_const - ELSE - phi = ATAN(my/mx) - pi_const - ENDIF - ENDIF - - ! write(36,'(i4,2f12.6)') mod(imesh,33),rho_11,rho_22 - ris(imesh,1) = rho_up - ris(imesh,2) = rho_down - den%theta_pw(imesh) = theta - den%phi_pw(imesh) = phi - ENDDO - - DO jspin = 1,input%jspins - fftwork=0.0 - CALL fft3d(ris(:,jspin),fftwork,den%pw(:,jspin),stars,-1) - ENDDO - - IF (.NOT.input%film) RETURN - - - !Now the vacuum part starts - - - ALLOCATE(rvacxy(ifft2,vacuum%nmzxyd,2,4)) - ALLOCATE (rz(vacuum%nmzd,2,2)) - !---> fouriertransform the diagonal part of the density matrix - !---> in the vacuum, rz & rxy, to real space (rvacxy) - DO iden = 1,2 - DO ivac = 1,vacuum%nvac - DO imz = 1,vacuum%nmzxyd - rziw = 0.0 - IF (oneD%odi%d1) THEN - CALL judft_error("oneD not implemented",calledby="rhodirgen") - !CALL fft2d(oneD%k3,odi%M,odi%n2d,rvacxy(0,imz,ivac,iden),fftwork,& - ! den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),& - ! vacuum,odi%nq2,odi%kimax2,1,& - ! & %igf,odl%pgf,odi%nst2) - ELSE - CALL fft2d(stars,rvacxy(:,imz,ivac,iden),fftwork,& - den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),& - vacuum%nmzxyd,1) - ENDIF - ENDDO - ENDDO - ENDDO - !---> fouriertransform the off-diagonal part of the density matrix - DO ivac = 1,vacuum%nvac - DO imz = 1,vacuum%nmzxyd - rziw = 0.0 - vz_r = den%vacz(imz,ivac,3) - vz_i = den%vacz(imz,ivac,4) - IF (oneD%odi%d1) THEN - CALL judft_error("oneD not implemented",calledby="rhodirgen") - !CALL fft2d(oneD%k3,odi%M,odi%n2d,& - ! rvacxy(0,imz,ivac,3),rvacxy(0,imz,ivac,4),& - ! vz_r,vz_i,den%vacxy(imz,1,ivac,3),& - ! vacuum,odi%nq2,odi%kimax2,1,& - ! & %igf,odl%pgf,odi%nst2) - ELSE - CALL fft2d(stars,rvacxy(:,imz,ivac,3),rvacxy(:,imz,ivac,4),& - vz_r,vz_i,den%vacxy(imz,1,ivac,3),vacuum%nmzxyd,1) - ENDIF - ENDDO - ENDDO - - !---> calculate the four components of the matrix potential on - !---> real space mesh - DO ivac = 1,vacuum%nvac - DO imz = 1,vacuum%nmzxyd - DO imesh = 1,ifft2 - rho_11 = rvacxy(imesh,imz,ivac,1) - rho_22 = rvacxy(imesh,imz,ivac,2) - rho_21r = rvacxy(imesh,imz,ivac,3) - rho_21i = rvacxy(imesh,imz,ivac,4) - mx = 2*rho_21r - my = -2*rho_21i - mz = (rho_11-rho_22) - magmom = SQRT(mx**2 + my**2 + mz**2) - rhotot = rho_11 + rho_22 - rho_up = (rhotot + magmom)/2 - rho_down= (rhotot - magmom)/2 - - IF (ABS(mz) .LE. eps) THEN - theta = pi_const/2 - ELSEIF (mz .GE. 0.0) THEN - theta = ATAN(SQRT(mx**2 + my**2)/mz) - ELSE - theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const - ENDIF - - IF (ABS(mx) .LE. eps) THEN - IF (ABS(my) .LE. eps) THEN - phi = 0.0 - ELSEIF (my .GE. 0.0) THEN - phi = pi_const/2 - ELSE - phi = -pi_const/2 - ENDIF - ELSEIF (mx .GE. 0.0) THEN - phi = ATAN(my/mx) - ELSE - IF (my .GE. 0.0) THEN - phi = ATAN(my/mx) + pi_const - ELSE - phi = ATAN(my/mx) - pi_const - ENDIF - ENDIF - - rvacxy(imesh,imz,ivac,1) = rho_up - rvacxy(imesh,imz,ivac,2) = rho_down - den%theta_vacxy(imesh,imz,ivac) = theta - den%phi_vacxy(imesh,imz,ivac) = phi - ENDDO - ENDDO - DO imz = vacuum%nmzxyd+1,vacuum%nmzd - rho_11 = den%vacz(imz,ivac,1) - rho_22 = den%vacz(imz,ivac,2) - rho_21r = den%vacz(imz,ivac,3) - rho_21i = den%vacz(imz,ivac,4) - mx = 2*rho_21r - my = -2*rho_21i - mz = (rho_11-rho_22) - magmom = SQRT(mx**2 + my**2 + mz**2) - rhotot = rho_11 + rho_22 - rho_up = (rhotot + magmom)/2 - rho_down= (rhotot - magmom)/2 - - IF (ABS(mz) .LE. eps) THEN - theta = pi_const/2 - ELSEIF (mz .GE. 0.0) THEN - theta = ATAN(SQRT(mx**2 + my**2)/mz) - ELSE - theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const - ENDIF - - IF (ABS(mx) .LE. eps) THEN - IF (ABS(my) .LE. eps) THEN - phi = 0.0 - ELSEIF (my .GE. 0.0) THEN - phi = pi_const/2 - ELSE - phi = -pi_const/2 - ENDIF - ELSEIF (mx .GE. 0.0) THEN - phi = ATAN(my/mx) - ELSE - IF (my .GE. 0.0) THEN - phi = ATAN(my/mx) + pi_const - ELSE - phi = ATAN(my/mx) - pi_const - ENDIF - ENDIF - - den%vacz(imz,ivac,1) = rho_up - den%vacz(imz,ivac,2) = rho_down - den%theta_vacz(imz,ivac) = theta - den%phi_vacz(imz,ivac) = phi - ENDDO - ENDDO - !---> Fouriertransform the matrix potential back to reciprocal space - DO jspin = 1,input%jspins - DO ivac = 1,vacuum%nvac - DO imz = 1,vacuum%nmzxyd - fftwork=0.0 - IF (oneD%odi%d1) THEN - call judft_error("oneD not implemented",calledby="rhodirgen") - !CALL fft2d(oneD%k3,odi%M,odi%n2d,& - ! rvacxy(0,imz,ivac,jspin),fftwork,& - ! den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),& - ! vacuum,odi%nq2,odi%kimax2,-1,& - ! & %igf,odl%pgf,odi%nst2) - ELSE - CALL fft2d(stars,rvacxy(:,imz,ivac,jspin),fftwork,& - den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),& - vacuum%nmzxyd,-1) - END IF - ENDDO - ENDDO - ENDDO - - RETURN - END SUBROUTINE rhodirgen -END MODULE m_rhodirgen diff --git a/vgen/visxcg.f90 b/vgen/visxcg.f90 deleted file mode 100644 index dff91b53d9764a19a49642d4ae19e8ab9dc7593a..0000000000000000000000000000000000000000 --- a/vgen/visxcg.f90 +++ /dev/null @@ -1,433 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- -MODULE m_visxcg - USE m_juDFT - ! ****************************************************** - ! subroutine generates the exchange-correlation potential - ! in the interstitial region c.l.fu - ! including gradient corrections. t.a. 1996. - ! ****************************************************** -CONTAINS - SUBROUTINE visxcg(ifftd,stars,sym, ifftxc3d, cell,den,& - xcpot,input, obsolete,noco, vxc,vx,exc) - - ! ****************************************************** - ! instead of visxcor.f: the different exchange-correlation - ! potentials defined through the key icorr are called through - ! the driver subroutine vxcallg.f,for the energy density - excallg - ! subroutines vectorized - ! ** r.pentcheva 22.01.96 - ! ********************************************************* - ! in case of total = .true. calculates the ex-corr. energy - ! density - ! ** r.pentcheva 08.05.96 - ! ****************************************************************** - - USE m_grdrsis - USE m_prpxcfftmap - USE m_mkgxyz3 - USE m_fft3d - USE m_fft3dxc - USE m_types - IMPLICIT NONE - - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_obsolete),INTENT(IN) :: obsolete - TYPE(t_input),INTENT(IN) :: input - TYPE(t_noco),INTENT(IN) :: noco - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_stars),INTENT(IN) :: stars - TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_potden),INTENT(IN) :: den - TYPE(t_potden),INTENT(INOUT) :: vxc,vx,exc - - ! .. - ! .. Scalar Arguments .. - ! .. - ! .. Array Arguments .. - - REAL rhmni ,d_15,sprsv - ! - !-----> fft information for xc potential + energy - ! - INTEGER, ALLOCATABLE :: igxc_fft(:) - REAL, ALLOCATABLE :: gxc_fft(:,:) - ! .. - ! .. Local Scalars .. - INTEGER :: i ,k,js,nt,ifftxc3,idm,jdm,ndm,ig - COMPLEX :: ci - REAL :: rhotot - INTEGER :: ifftd,ifftxc3d - TYPE(t_gradients)::grad - - ! .. - ! .. Local Arrays .. - COMPLEX, ALLOCATABLE :: fg3(:),cqpw(:,:), ph_wrk(:) - REAL, ALLOCATABLE :: bf3(:) - REAL, ALLOCATABLE :: rho(:,:),rhd1(:,:,:),rhd2(:,:,:) - REAL, ALLOCATABLE :: mx(:),my(:) - REAL, ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:) - ! - REAL, ALLOCATABLE :: v_x(:,:),v_xc(:,:),e_xc(:),vcon(:) - - ! .. unused input (needed for other noco GGA-implementations) .. - - !ta+ - !.....------------------------------------------------------------------ - !-------> abbreviations - ! - ! ph_wrk: work array containing phase * g_x,gy...... - ! den%pw: charge density stored as stars - ! rho : charge density stored in real space - ! v_xc : exchange-correlation potential in real space - ! exc : exchange-correlation energy density in real space - ! kxc1d : dimension of the charge density fft box in the pos. domain - ! kxc2d : defined in dimens.f program (subroutine apws).1,2,3 indic - ! kxc3d ; a_1, a_2, a_3 directions. - ! kq(i) : i=1,2,3 actual length of the fft-box for which fft is done - ! nstr : number of members (arms) of reciprocal lattice (g) vector - ! of each star. - ! nxc3_fft: number of stars in the charge density fft-box - ! ng3 : number of 3 dim. stars in the charge density sphere define - ! by gmax - ! kmxxc_fft: number of g-vectors forming the nxc3_fft stars in the - ! charge density or xc-density sphere - ! kimax : number of g-vectors forming the ng3 stars in the gmax-sphe - ! ifftxc3d: elements (g-vectors) in the charge density fft-box - ! igfft : pointer from the g-sphere (stored as stars) to fft-grid - ! and from fft-grid to g-sphere (stored as stars) - ! pgfft : contains the phases of the g-vectors of sph. - ! isn : isn = +1, fft transform for g-space to r-space - ! isn = -1, vice versa - ! - !------------------------------------------------------------------- - ! - !---> set up pointer for backtransformation of from g-vector in - ! positive domain of xc density fftbox into stars. - ! also the x,y,z components of the g-vectors are set up to calculate - ! derivatives. - ! in principle this can also be done in main program once. - ! it is done here to save memory. - ! - ifftd=27*stars%mx1*stars%mx2*stars%mx3 - ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft - ALLOCATE ( igxc_fft(0:ifftxc3d-1),gxc_fft(0:ifftxc3d-1,3) ) - CALL prp_xcfft_map(stars,sym, cell, igxc_fft,gxc_fft) - ! - ifftxc3=stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft - - IF (stars%ng3.GT.stars%ng3) THEN - WRITE(6,'(/'' stars%ng3.gt.stars%ng3. stars%ng3,stars%ng3='',2i6)') stars%ng3,stars%ng3 - CALL juDFT_error("ng3.gt.n3d",calledby="visxcg") - ENDIF - - d_15=1.e-15 - ! - ci=CMPLX(0.,1.) - ! - ! Allocate arrays - ! ff - ALLOCATE( bf3(0:ifftd-1),ph_wrk(0:ifftxc3d-1), rho(0:ifftxc3d-1,input%jspins),& - rhd1(0:ifftxc3d-1,input%jspins,3), rhd2(0:ifftxc3d-1,input%jspins,6) ) - IF (noco%l_noco) THEN - ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),& - magmom(0:ifftxc3-1), dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) ) - END IF - - - !--> transform charge density to real space - - DO js=1,input%jspins - CALL fft3dxc(rho(0:,js),bf3, den%pw(:,js), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) - END DO - - IF (noco%l_noco) THEN - - ! for off-diagonal parts the same - CALL fft3dxc(mx,my, den%pw(:,3), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) - - DO i=0,ifftxc3-1 - rhotot= 0.5*( rho(i,1) + rho(i,2) ) - magmom(i)= SQRT( (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 ) - rho(i,1)= rhotot+magmom(i) - rho(i,2)= rhotot-magmom(i) - END DO - - ENDIF - - IF (.not.xcpot%is_gga()) GOTO 100 - - ! In collinear calculations all derivatives are calculated in g-spce, - ! in non-collinear calculations the derivatives of |m| are calculated in real space. - - !--> for d(rho)/d(x,y,z) = rhd1(:,:,idm) (idm=1,2,3). - ! - ! ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z). - - ALLOCATE(cqpw(stars%ng3,input%jspins)) - - DO js= 1,input%jspins - DO i = 1,stars%ng3 - cqpw(i,js)= ci*den%pw(i,js) - END DO - END DO - - DO idm=1,3 - - DO ig = 0 , stars%kmxxc_fft - 1 - ph_wrk(ig) = stars%pgfft(ig) * gxc_fft(ig,idm) - END DO - - DO js=1,input%jspins - CALL fft3dxc(rhd1(0:,js,idm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,& - stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr) - END DO - - END DO - - IF (noco%l_noco) THEN - - CALL grdrsis(magmom,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete, dmagmom ) - - DO i=0,ifftxc3-1 - DO idm=1,3 - rhotot= rhd1(i,1,idm)/2.+rhd1(i,2,idm)/2. - rhd1(i,1,idm)= rhotot+dmagmom(i,idm) - rhd1(i,2,idm)= rhotot-dmagmom(i,idm) - END DO - END DO - - END IF - - - !--> for dd(rho)/d(xx,xy,yy,zx,yz,zz) = rhd2(:,:,idm) (idm=1,2,3,4,5,6) - ! - ! ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z) * g_(x,y,z) - - DO i = 1,stars%ng3 - DO js=1,input%jspins - cqpw(i,js)= -den%pw(i,js) - END DO - END DO - - ndm = 0 - DO idm = 1,3 - DO jdm = 1,idm - ndm = ndm + 1 - - DO ig = 0 , stars%kmxxc_fft-1 - ph_wrk(ig) = stars%pgfft(ig)*gxc_fft(ig,idm)*gxc_fft(ig,jdm) - ENDDO - - DO js=1,input%jspins - CALL fft3dxc(rhd2(0:,js,ndm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,& - stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr) - END DO - END DO ! jdm - END DO ! idm - - DEALLOCATE(cqpw) - - IF (noco%l_noco) THEN - - DO idm = 1,3 - CALL grdrsis(dmagmom(0,idm),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete, ddmagmom(0,1,idm) ) - END DO - - ndm= 0 - DO idm = 1,3 - DO jdm = 1,idm - ndm = ndm + 1 - - DO i=0,ifftxc3-1 - rhotot= rhd2(i,1,ndm)/2.+rhd2(i,2,ndm)/2. - rhd2(i,1,ndm)= rhotot + ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. - rhd2(i,2,ndm)= rhotot - ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. - END DO - - ENDDO !jdm - ENDDO !idm - - END IF - -100 CONTINUE - - - DEALLOCATE ( ph_wrk ) - IF (noco%l_noco) THEN - DEALLOCATE(mx,my,magmom,dmagmom,ddmagmom) - END IF - ! - DO js=1,input%jspins - DO i=0,ifftxc3-1 - rho(i,js)=MAX(rho(i,js),d_15) - ENDDO - END DO - bf3=0.0 - ! allocate the other arrays - ! - CALL xcpot%alloc_gradients(ifftxc3d,input%jspins,grad) - - ! - ! calculate the quantities such as abs(grad(rho)),.. used in - ! evaluating the gradient contributions to potential and energy. - ! - CALL mkgxyz3 (ifftxc3d,input%jspins,ifftxc3,input%jspins,rho, rhd1(0,1,1),rhd1(0,1,2),rhd1(0,1,3),& - rhd2(0,1,1),rhd2(0,1,3),rhd2(0,1,6), rhd2(0,1,5),rhd2(0,1,4),rhd2(0,1,2), grad) - - DEALLOCATE ( rhd1,rhd2 ) - ALLOCATE ( v_xc(0:ifftxc3d-1,input%jspins) ) - ALLOCATE ( v_x (0:ifftxc3d-1,input%jspins) ) - ! - ! calculate the exchange-correlation potential in real space - ! - nt=ifftxc3 - ! - ! rhmni: rho_minimum_interstitial. - - rhmni=10.e+10 - - DO js=1,input%jspins - DO i=0,ifftxc3-1 - rho(i,js)=MAX(rho(i,js),d_15) - rhmni=MIN(rhmni,rho(i,js)) - ENDDO - ENDDO - - - IF (rhmni.LT.obsolete%chng) THEN - WRITE(6,'(/'' rhmn.lt.obsolete%chng in visxc. rhmn,obsolete%chng='',2d9.2)') rhmni,obsolete%chng - ! CALL juDFT_error("visxcg: rhmn.lt.chng",calledby="visxcg") - ENDIF - CALL xcpot%get_vxc(input%jspins,rho,v_xc,v_x,grad) - ! - !----> back fft to g space - !----> perform back fft transform: v_xc(r) --> vxc(star) - ! - ALLOCATE(fg3(stars%ng3)) - - DO js = 1,input%jspins - bf3=0.0 - CALL fft3dxc(v_xc(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) - ! - DO k = 1,stars%nxc3_fft - vxc%pw(k,js) = vxc%pw(k,js) + fg3(k) - ENDDO - - ! - !====> INCLUDING TOTAL ENERGY - ! - IF (input%total) THEN - ! - !----> Perform fft transform: v_xc(star) --> vxc(r) - ! !Use large fft mesh for convolution - ! - fg3(stars%nxc3_fft+1:)=0.0 - ALLOCATE ( vcon(0:ifftd-1) ) - CALL fft3d(vcon(0),bf3, fg3, stars,+1) - ! - !----> Convolute with step function - ! - DO i=0,ifftd-1 - vcon(i)=stars%ufft(i)*vcon(i) - ENDDO - bf3=0.0 - CALL fft3d(vcon(0),bf3, fg3, stars,-1,.FALSE.) - DEALLOCATE ( vcon ) - ! - !----> add to warped coulomb potential - ! - DO k = 1,stars%ng3 - vxc%pw_w(k,js) = vxc%pw_w(k,js) + fg3(k) - ENDDO - - ENDIF - - bf3=0.0 - CALL fft3dxc(v_x(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft, & - stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) - ! - DO k = 1,stars%nxc3_fft - vx%pw(k,js) = vx%pw(k,js) + fg3(k) - ENDDO - - ! - !====> INCLUDING TOTAL ENERGY - ! - IF (input%total) THEN - ! - !----> Perform fft transform: v_xc(star) --> vxc(r) - ! !Use large fft mesh for convolution - ! - - fg3(stars%nxc3_fft+1:)=0.0 - ALLOCATE ( vcon(0:ifftd-1) ) - CALL fft3d(vcon(0),bf3, fg3, stars,+1) - ! - !----> Convolute with step function - ! - DO i=0,ifftd-1 - vcon(i)=stars%ufft(i)*vcon(i) - ENDDO - bf3=0.0 - CALL fft3d(vcon(0),bf3, fg3, stars,-1,.FALSE.) - DEALLOCATE ( vcon ) - ! - !----> add to warped exchange-potential - ! - DO k = 1,stars%ng3 - vx%pw_w(k,js) = vx%pw_w(k,js) + fg3(k) - ENDDO - - ENDIF - - ENDDO - DEALLOCATE ( v_x,v_xc ) - ! - ! calculate the ex.-cor energy density in real space - ! - IF (ALLOCATED(exc%pw_w)) THEN - ALLOCATE ( e_xc(0:ifftxc3d-1) ) - CALL xcpot%get_exc(input%jspins,rho,e_xc,grad) - ! - !----> perform back fft transform: exc(r) --> exc(star) - ! - bf3=0.0 - CALL fft3dxc(e_xc,bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) - DEALLOCATE ( e_xc ) - ! - !----> Perform fft transform: exc(star) --> exc(r) - ! !Use large fft mesh for convolution - ! - fg3(stars%nxc3_fft+1:)=0.0 - bf3=0.0 - ALLOCATE ( vcon(0:ifftd-1) ) - CALL fft3d(vcon,bf3,fg3, stars,+1) - - DO i=0,ifftd-1 - vcon(i)=stars%ufft(i)*vcon(i) - ENDDO - ! - ! ---> back fft to g space - ! - bf3=0.0 - CALL fft3d(vcon,bf3,exc%pw_w(:,1), stars,-1,.FALSE.) - DEALLOCATE ( vcon ) - ! - ENDIF - - DEALLOCATE(fg3) - DEALLOCATE ( bf3,rho,igxc_fft,gxc_fft ) - - - - - END SUBROUTINE visxcg -END MODULE m_visxcg diff --git a/vgen/vmtxc.f90 b/vgen/vmtxc.f90 deleted file mode 100644 index 77e0e081e392916020223fe20fb5c3ab5d497bc4..0000000000000000000000000000000000000000 --- a/vgen/vmtxc.f90 +++ /dev/null @@ -1,145 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- -MODULE m_vmtxc -CONTAINS - SUBROUTINE vmtxc( DIMENSION,sphhar,atoms,den,xcpot,input,sym, vxc, exc,vx) - ! ******************************************************************** - ! Calculate the LDA xc-potential in the MT-spheres - ! ********************************************************* - - USE m_lhglpts - USE m_gaussp - USE m_types - IMPLICIT NONE - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_dimension),INTENT(IN) :: DIMENSION - TYPE(t_input),INTENT(IN) :: input - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_sphhar),INTENT(IN) :: sphhar - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_potden),INTENT(IN) :: den - TYPE(t_potden),INTENT(INOUT) :: vxc,exc - TYPE(t_potden),INTENT(INOUT),optional :: vx - - ! .. - ! .. Local Scalars .. - REAL elh,rr2,r2rho,vlh - INTEGER jr,js,k,lh,n,nd,i,nsp,nat - ! .. - ! .. Local Arrays .. - REAL v_x(DIMENSION%nspd,DIMENSION%jspd),v_xc(DIMENSION%nspd,DIMENSION%jspd),pos0(3) - REAL rhoxc(DIMENSION%nspd,DIMENSION%jspd),e_ex(DIMENSION%nspd),wt(DIMENSION%nspd),rx(3,DIMENSION%nspd) - REAL ylh(DIMENSION%nspd,0:sphhar%nlhd,sphhar%ntypsd) - ! .. - ! generates nspd points on a sherical shell with radius 1.0 - ! angular mesh equidistant in phi, - ! theta are zeros of the legendre polynomials - ! - CALL gaussp(atoms%lmaxd,rx,wt) - - nsp = DIMENSION%nspd - ! - ! generates the lattice harmonics on the angular mesh - ! - CALL lhglpts( sphhar,atoms, rx,nsp, sym, ylh) - ! - ! - ! loop over topologically non-equivalent atoms - ! - !$OMP PARALLEL DO DEFAULT(NONE)& - !$OMP SHARED(atoms,input,sphhar,xcpot,den,ylh,nsp,wt,vxc,vx,exc)& - !$OMP PRIVATE(elh,rr2,r2rho,vlh,jr,js,k,lh,n,nd,i,nat)& - !$OMP PRIVATE(v_x,v_xc,pos0,rhoxc,e_ex) - DO n = 1,atoms%ntype - nat=SUM(atoms%neq(:n-1))+1 - nd = atoms%ntypsy(nat) - ! - ! loop over radial mesh - ! - DO jr = 1,atoms%jri(n) - rr2 = 1.e0/ (atoms%rmsh(jr,n)*atoms%rmsh(jr,n)) - DO js = 1,input%jspins - rhoxc(:,js) = 0. - DO lh = 0,sphhar%nlh(nd) - r2rho = rr2*den%mt(jr,lh,n,js) - ! - ! generate the densities on an angular mesh - ! - DO k = 1,nsp - rhoxc(k,js) = rhoxc(k,js) + ylh(k,lh,nd)*r2rho - ENDDO - ENDDO - ENDDO - ! calculate the ex.-cor. potential - ! - CALL xcpot%get_vxc(input%jspins,rhoxc(:nsp,:), v_xc,v_x) - ! ----> now determine the corresponding potential number - DO js = 1,input%jspins - ! - ! ---> multiplikate vxc with the weights of the k-points - ! - DO k = 1,nsp - v_xc(k,js) = v_xc(k,js)*wt(k) - - v_x (k,js) = v_x (k,js)*wt(k) - ENDDO - DO lh = 0,sphhar%nlh(nd) - - ! - ! ---> determine the corresponding potential number through gauss integration - ! - vlh=DOT_PRODUCT( v_xc(:nsp,js),ylh(:nsp,lh,nd)) - - ! - ! ---> add to the given potential - - vxc%mt(jr,lh,n,js) = vxc%mt(jr,lh,n,js) + vlh - - IF (PRESENT(vx) ) THEN - vlh = 0 - DO k = 1,nsp - vlh = vlh + v_x(k,js)*ylh(k,lh,nd) - END DO - - vx%mt(jr,lh,n,js) = vx%mt(jr,lh,n,js) + vlh - ENDIF - ENDDO - ENDDO - IF (ALLOCATED(exc%mt)) THEN - ! - ! calculate the ex.-cor energy density - ! - CALL xcpot%get_exc(input%jspins,rhoxc(:nsp,:), e_ex) - ! ----> now determine the corresponding energy density number - ! - ! ---> multiplikate exc with the weights of the k-points - ! - DO k = 1,nsp - e_ex(k) = e_ex(k)*wt(k) - ENDDO - DO lh = 0,sphhar%nlh(nd) - - ! - ! ---> determine the corresponding potential number through gauss integration - ! - elh=DOT_PRODUCT(e_ex(:nsp),ylh(:nsp,lh,nd)) - - ! - ! ---> add to the given potential - - exc%mt(jr,lh,n,1) = elh - ENDDO - END IF - - ENDDO - - ENDDO - !$OMP END PARALLEL DO - ! - - RETURN - END SUBROUTINE vmtxc -END MODULE m_vmtxc diff --git a/vgen/vmtxcg.F90 b/vgen/vmtxcg.F90 deleted file mode 100644 index 00416a376d375afb4129ba9cb1a4e236f9bd58d9..0000000000000000000000000000000000000000 --- a/vgen/vmtxcg.F90 +++ /dev/null @@ -1,381 +0,0 @@ -!-------------------------------------------------------------------------------- -! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany -! This file is part of FLEUR and available as free software under the conditions -! of the MIT license as expressed in the LICENSE file in more detail. -!-------------------------------------------------------------------------------- -MODULE m_vmtxcg - !.....------------------------------------------------------------------ - ! Calculate the GGA xc-potential in the MT-spheres - !.....------------------------------------------------------------------ - ! instead of vmtxcor.f: the different exchange-correlation - ! potentials defined through the key icorr are called through - ! the driver subroutine vxcallg.f, subroutines vectorized - ! ** r.pentcheva 22.01.96 - ! ********************************************************* - ! angular mesh calculated on speacial gauss-legendre points - ! in order to use orthogonality of lattice harmonics and - ! avoid a least square fit - ! ** r.pentcheva 04.03.96 - ! ********************************************************* - ! MPI and OpenMP parallelization - ! U.Alekseeva, February 2017 - ! ********************************************************* - -CONTAINS - SUBROUTINE vmtxcg(dimension,mpi,sphhar,atoms,& - den,xcpot,input,sym, obsolete,vxc,vx,exc) - -#include"cpp_double.h" - - USE m_grdchlh - USE m_mkgylm - USE m_gaussp - USE m_types_xcpot_inbuild - USE m_types - IMPLICIT NONE - - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_dimension),INTENT(IN) :: dimension - TYPE(t_mpi),INTENT(IN) :: mpi - TYPE(t_obsolete),INTENT(IN) :: obsolete - TYPE(t_input),INTENT(IN) :: input - TYPE(t_sym),INTENT(IN) :: sym - TYPE(t_sphhar),INTENT(IN) :: sphhar - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_potden),INTENT(IN) :: den - TYPE(t_potden),INTENT(INOUT) :: vxc,vx,exc -#ifdef CPP_MPI - include "mpif.h" -#endif - ! .. - ! .. Local Scalars .. - TYPE(t_gradients)::grad - INTEGER jr,js,k,lh,n,nd,ist,nsp,ixpm ,i,nat - REAL rhmnm,d_15,elh,vlh - LOGICAL lwbc ! if true, white-bird trick - ! .. - ! .. Local Arrays .. - REAL v_x(dimension%nspd,dimension%jspd),v_xc(dimension%nspd,dimension%jspd),e_xc(dimension%nspd),rx(3,dimension%nspd) - REAL vxcl(DIMENSION%nspd,DIMENSION%jspd),excl(DIMENSION%nspd),divi - TYPE(t_xcpot_inbuild)::xcpot_tmp - REAL wt(dimension%nspd),rr2(atoms%jmtd),thet(dimension%nspd) - REAL, ALLOCATABLE :: ylh(:,:,:),ylht(:,:,:),ylhtt(:,:,:) - REAL, ALLOCATABLE :: ylhf(:,:,:),ylhff(:,:,:),ylhtf(:,:,:) - REAL, ALLOCATABLE :: chlh(:,:,:),chlhdr(:,:,:),chlhdrr(:,:,:) - REAL, ALLOCATABLE :: ch(:,:),chdr(:,:),chdt(:,:),chdf(:,:) - REAL, ALLOCATABLE :: chdrr(:,:),chdtt(:,:),chdff(:,:),chdtf(:,:) - REAL, ALLOCATABLE :: chdrt(:,:),chdrf(:,:) - - !locals for mpi - integer :: ierr - integer:: n_start,n_stride -#ifdef CPP_MPI - REAL :: vr_local(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) - REAL :: vxr_local(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) - REAL :: excr_local(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) -#endif - ! .. - ! .. Intrinsic Functions .. - INTRINSIC max,mod,min - !ta+ - !.....------------------------------------------------------------------ - -#ifdef CPP_MPI - CALL MPI_BCAST(obsolete%ndvgrd,1,MPI_INTEGER,0,mpi%mpi_comm,ierr) - vr_local = 0.d0 - vxr_local = 0.d0 - excr_local = 0.d0 -#endif - - lwbc=.false. - - d_15 = 1.e-15 - ixpm = 1 - ist = 1 - - ! generates nspd points on a sherical shell with radius 1.0 - ! angular mesh equidistant in phi, - ! theta are zeros of the legendre polynomials - ! - CALL gaussp(atoms%lmaxd, rx,wt) - - nsp = dimension%nspd - ! - ! generates the lattice harmonics on the angular mesh - ! - ALLOCATE ( ylh(dimension%nspd,0:sphhar%nlhd,sphhar%ntypsd),ylht(dimension%nspd,0:sphhar%nlhd,sphhar%ntypsd),& - & ylhtt(dimension%nspd,0:sphhar%nlhd,sphhar%ntypsd),ylhf(dimension%nspd,0:sphhar%nlhd,sphhar%ntypsd),& - & ylhff(dimension%nspd,0:sphhar%nlhd,sphhar%ntypsd),ylhtf(dimension%nspd,0:sphhar%nlhd,sphhar%ntypsd) ) - - CALL lhglptg(& - & sphhar,atoms,& - & rx,nsp,xcpot,sym,& - & ylh,thet,ylht,ylhtt,ylhf,ylhff,ylhtf) - ! - ! - ! loop over topologically non-equivalent atoms - ! - - -#ifdef CPP_MPI - n_start=mpi%irank+1 - n_stride=mpi%isize -#else - n_start=1 - n_stride=1 -#endif - - ALLOCATE ( chlh(atoms%jmtd,0:sphhar%nlhd,dimension%jspd),& - & chlhdr(atoms%jmtd,0:sphhar%nlhd,dimension%jspd),chlhdrr(atoms%jmtd,0:sphhar%nlhd,dimension%jspd)) - - DO 200 n = n_start,atoms%ntype,n_stride - SELECT TYPE(xcpot) - TYPE IS(t_xcpot_inbuild) - IF (xcpot%lda_atom(n))THEN - IF((.NOT.xcpot%is_name("pw91"))) & - CALL judft_warn("Using locally LDA only possible with pw91 functional") - CALL xcpot_tmp%init("l91",.FALSE.,atoms%ntype) - ENDIF - END SELECT - - nat=sum(atoms%neq(:n-1))+1 - nd = atoms%ntypsy(nat) - - DO jr = 1,atoms%jri(n) - rr2(jr) = 1.e0/ (atoms%rmsh(jr,n)*atoms%rmsh(jr,n)) - ENDDO - - DO lh = 0,sphhar%nlh(nd) - - ! calculates gradients of radial charge densities of l=> 0. - ! rho*ylh/r**2 is charge density. chlh=rho/r**2. - ! charge density=sum(chlh*ylh). - ! chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr. - - DO js = 1,input%jspins - - DO jr = 1,atoms%jri(n) - chlh(jr,lh,js) = den%mt(jr,lh,n,js)*rr2(jr) - ENDDO - - CALL grdchlh(ixpm,ist,atoms%jri(n),atoms%dx(n),atoms%rmsh(1,n),& - chlh(1,lh,js),obsolete%ndvgrd, chlhdr(1,lh,js),chlhdrr(1,lh,js)) - - ENDDO ! js - ENDDO ! lh - - ! - !--> loop over radial mesh - ! - !!$OMP PARALLEL DEFAULT(none) & -#ifdef CPP_MPI - !!$OMP& SHARED(vr_local,vxr_local,excr_local) & -#endif - !!$OMP& SHARED(vxc,vx,exc) & - !!$OMP& SHARED(dimension,mpi,sphhar,atoms,den,xcpot,input,sym,obsolete)& - !!$OMP& SHARED(n,nd,ist,ixpm,nsp,nat,d_15,lwbc) & - !!$OMP& SHARED(rx,wt,rr2,thet,xcpot_tmp) & - !!$OMP& SHARED(ylh,ylht,ylhtt,ylhf,ylhff,ylhtf) & - !!$OMP& SHARED(chlh,chlhdr,chlhdrr) & - !!$OMP& SHARED(ierr,n_start,n_stride) & - !!$OMP& PRIVATE(js,k,lh,i,rhmnm,elh,vlh) & - !!$OMP& PRIVATE(v_x,v_xc,e_xc,excl,vxcl,divi) & - !!$OMP& PRIVATE(grad) & - !!$OMP& PRIVATE(ch,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf) - ALLOCATE ( ch(DIMENSION%nspd,DIMENSION%jspd),chdr(DIMENSION%nspd,DIMENSION%jspd),& - chdt(DIMENSION%nspd,DIMENSION%jspd), chdf(DIMENSION%nspd,DIMENSION%jspd),& - chdrr(DIMENSION%nspd,DIMENSION%jspd),chdtt(DIMENSION%nspd,DIMENSION%jspd), & - chdff(DIMENSION%nspd,DIMENSION%jspd),chdtf(DIMENSION%nspd,DIMENSION%jspd),& - chdrt(DIMENSION%nspd,DIMENSION%jspd), chdrf(DIMENSION%nspd,DIMENSION%jspd)) - !$OMP DO - DO 190 jr = 1,atoms%jri(n) - ! - ! following are at points on jr-th sphere. - ! initialize arrays - ! - ch(:,:) = 0.0 ! charge density - chdr(:,:) = 0.0 ! d(ch)/dr - chdt(:,:) = 0.0 ! d(ch)/dtheta - chdf(:,:) = 0.0 ! d(ch)/dfai - chdrr(:,:) = 0.0 ! dd(ch)/drr - chdtt(:,:) = 0.0 ! dd(ch)/dtt - chdff(:,:) = 0.0 ! dd(ch)/dff - chdtf(:,:) = 0.0 ! dd(ch)/dtf - chdrt(:,:) = 0.0 ! d(d(ch)/dr)dt - chdrf(:,:) = 0.0 ! d(d(ch)/dr)df - - DO js = 1,input%jspins - DO lh = 0,sphhar%nlh(nd) - - ! generate the densities on an angular mesh - - DO k = 1,nsp - ch(k,js) = ch(k,js) + ylh(k,lh,nd)*chlh(jr,lh,js) - ENDDO - - ! - DO k = 1,nsp - chdr(k,js) =chdr(k,js)+ ylh(k,lh,nd)*chlhdr(jr,lh,js) - chdrr(k,js)=chdrr(k,js)+ylh(k,lh,nd)*chlhdrr(jr,lh,js) - ENDDO - - DO k = 1,nsp - chdrt(k,js)=chdrt(k,js)+ylht(k,lh,nd)*chlhdr(jr,lh,js) - chdrf(k,js)=chdrf(k,js)+ylhf(k,lh,nd)*chlhdr(jr,lh,js) - chdt(k,js) =chdt(k,js) +ylht(k,lh,nd)*chlh(jr,lh,js) - chdf(k,js) =chdf(k,js) +ylhf(k,lh,nd)*chlh(jr,lh,js) - chdtt(k,js)=chdtt(k,js)+ylhtt(k,lh,nd)*chlh(jr,lh,js) - chdff(k,js)=chdff(k,js)+ylhff(k,lh,nd)*chlh(jr,lh,js) - chdtf(k,js)=chdtf(k,js)+ylhtf(k,lh,nd)*chlh(jr,lh,js) - ENDDO - - - ENDDO ! lh - ENDDO ! js - - CALL xcpot%alloc_gradients(nsp,input%jspins,grad) - CALL mkgylm(input%jspins,atoms%rmsh(jr,n),thet,nsp,DIMENSION%nspd,DIMENSION%jspd,ch,chdr,& - chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf,grad) - - !Set charge to minimum value - ch(:,:)=MAX(ch(:,:),d_15) - ! - ! calculate the ex.-cor. potential - - CALL xcpot%get_vxc(input%jspins,ch(:nsp,:),v_xc,v_x,grad) - SELECT TYPE(xcpot) - TYPE IS (t_xcpot_inbuild) - IF (xcpot%lda_atom(n)) THEN - ! Use local part of pw91 for this atom - CALL xcpot_tmp%get_vxc(input%jspins,ch(:nsp,:),vxcl,v_x,grad) - !Mix the potentials - divi = 1.0 / (atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(1,n)) - v_xc(:,:) = ( vxcl(:,:) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +& - v_xc(:,:) * ( atoms%rmsh(jr,n) - atoms%rmsh(1,n) ) ) * divi - ENDIF - END SELECT - - IF (mpi%irank == 0) THEN - IF (mod(jr,1000).eq.0) WRITE (6,'(/'' 999vxc''/(10d15.7))')& - ((v_xc(k,js),k=1,nsp),js=1,input%jspins) - ENDIF !irank==0 - - - ! now determine the corresponding potential number - - DO js = 1,input%jspins - ! - ! multiplicate vx/vxc with the weights of the k-points - ! - DO k = 1,nsp - v_x (k,js) = v_x (k,js)*wt(k) - v_xc(k,js) = v_xc(k,js)*wt(k) - ENDDO - - IF (mpi%irank == 0) THEN - IF (mod(jr,1500).EQ.0) WRITE (6,'('' 999wt''/(10d15.7))') (wt(k),k=1,nsp) - IF (mod(jr,1500).EQ.0) WRITE (6,'('' 999vxc''/(10d15.7))') (v_xc(k,js),k=1,nsp) - ENDIF !irank==0 - - DO lh = 0,sphhar%nlh(nd) - ! - ! ---> determine the corresponding potential number - !c through gauss integration - ! - vlh=dot_product(v_xc(:nsp,js),ylh(:nsp,lh,nd)) -#ifdef CPP_MPI - vr_local(jr,lh,n,js) = vr_local(jr,lh,n,js) + vlh -#else - vxc%mt(jr,lh,n,js) = vxc%mt(jr,lh,n,js) + vlh -#endif - - ! ---> add to the given potential - - IF (mod(jr,1500).EQ.0) THEN - WRITE(6,'('' 999lh,js,jr,n,vlh='',4i5,d15.7)')& - & lh,js,jr,n,vlh - WRITE(6,'('' 9ylh''/(10d15.7))') (ylh(k,lh,nd),k=1,nsp) - ENDIF - - vlh=dot_product(v_x(:nsp,js),ylh(:nsp,lh,nd)) - -#ifdef CPP_MPI - vxr_local(jr,lh,n,js) = vxr_local(jr,lh,n,js) + vlh -#else - vx%mt(jr,lh,n,js) = vx%mt(jr,lh,n,js) + vlh -#endif - - ENDDO ! lh - ENDDO ! js - - IF (ALLOCATED(exc%mt)) THEN - ! - ! calculate the ex.-cor energy density - ! - CALL xcpot%get_exc(input%jspins,ch(:nsp,:),e_xc,grad) - SELECT TYPE(xcpot) - TYPE IS(t_xcpot_inbuild) - IF (xcpot%lda_atom(n)) THEN - ! Use local part of pw91 for this atom - CALL xcpot_tmp%get_exc(input%jspins,ch(:nsp,:),excl,grad) - !Mix the potentials - e_xc(:) = ( excl(:) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +& - e_xc(:) * ( atoms%rmsh(jr,n) - atoms%rmsh(1,n) ) ) * divi - ENDIF - END SELECT - IF (mpi%irank == 0) THEN - IF (mod(jr,10000).EQ.0) WRITE (6,'(/'' 999exc''/(10d15.7))') (e_xc(k),k=1,nsp) - ENDIF !irank==0 - - - ! now determine the corresponding energy density number - ! - ! multiplicate exc with the weights of the k-points - ! - DO k = 1,nsp - e_xc(k) = e_xc(k)*wt(k) - ENDDO - - DO lh = 0,sphhar%nlh(nd) - ! - ! determine the corresponding potential number through gauss - ! integration - ! - elh = 0.0 - DO k = 1,nsp - elh = elh + e_xc(k)*ylh(k,lh,nd) - ENDDO -#ifdef CPP_MPI - excr_local(jr,lh,n) = elh -#else - exc%mt(jr,lh,n,1) = elh -#endif - - ENDDO - ENDIF - -190 ENDDO - !!$OMP END DO - DEALLOCATE (ch,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf) - !!$OMP END PARALLEL - - ! WRITE(6,'(/'' n='',i3/'' 9999vr''/(10d15.7))') n, - ! & (((vr(jr,lh,n,js),jr=1,jri(n),100),lh=0,ntypsy(nat)),js=1,jspins) - ! WRITE(6,'(/'' 9999excr''/(10d15.7))') - ! & ((excr(jr,lh,n),jr=1,jri(n),100),lh=0,ntypsy(nat)) - -200 ENDDO - DEALLOCATE (chlh,chlhdr,chlhdrr) - -#ifdef CPP_MPI - CALL MPI_ALLREDUCE(vxr_local,vx%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,CPP_MPI_REAL,MPI_SUM,mpi%mpi_comm,ierr) !ToDo:CPP_MPI_REAL? - !using vxr_local as a temporal buffer - CALL MPI_ALLREDUCE(vr_local,vxr_local,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,CPP_MPI_REAL,MPI_SUM,mpi%mpi_comm,ierr) - vxc%mt = vxc%mt + vxr_local - CALL MPI_ALLREDUCE(excr_local,exc%mt(:,:,:,1),atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype,CPP_MPI_REAL,MPI_SUM,mpi%mpi_comm,ierr) -#endif - ! - - RETURN - END SUBROUTINE vmtxcg -END MODULE m_vmtxcg