diff --git a/PWCOND/do_cond.f90 b/PWCOND/do_cond.f90 index d161a4012..f248f6f3e 100644 --- a/PWCOND/do_cond.f90 +++ b/PWCOND/do_cond.f90 @@ -5,141 +5,136 @@ ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! -subroutine do_cond(nodenumber) +#include "f_defs.h" +! +SUBROUTINE do_cond(nodenumber) ! ! This is the main driver of the pwcond.x program. ! It calculates the complex band structure, solves the ! scattering problem and calculates the transmission coefficient. ! -#include "f_defs.h" USE ions_base, ONLY : nat, ityp, ntyp => nsp, tau - use pwcom - use cond - use io_files - use io_global, only : ionode_id -#ifdef __PARA - use para, only: me - use mp -#endif -implicit none - character(len=3) nodenumber - real(kind=DP), parameter :: eps=1.d-8 - real(kind=DP) :: wtot - integer :: ik, ien, ios, orbin, orbfin - logical :: write0 + USE pwcom + USE cond + USE io_files + USE io_global, ONLY : ionode, ionode_id - namelist /inputcond/ outdir, prefix, band_file, tran_file, fil_loc, & + USE mp + + IMPLICIT NONE + + CHARACTER(len=3) nodenumber + REAL(kind=DP), PARAMETER :: eps=1.d-8 + REAL(kind=DP) :: wtot + INTEGER :: ik, ien, ios, orbin, orbfin + LOGICAL :: write0 + + NAMELIST /inputcond/ outdir, prefix, band_file, tran_file, fil_loc, & lwrite_loc, lread_loc, ikind, iofspin, llocal, & bdl1, bdl2, bdr1, bdr2, nz1, dnslab, energy0, & denergy, ecut2d, ewind, epsproj, delgep,cutplot,& llapack - + CHARACTER (LEN=80) :: input_file - INTEGER :: nargs, iiarg, ierr, ilen + INTEGER :: nargs, iiarg, ierr, ILEN INTEGER, EXTERNAL :: iargc nd_nmbr=nodenumber - call init_clocks(.TRUE.) - call start_clock('PWCOND') - call start_clock('init') + CALL init_clocks(.TRUE.) + CALL start_clock('PWCOND') + CALL start_clock('init') ! ! set default values for variables in namelist ! - outdir = './' - prefix = ' ' - band_file = ' ' - tran_file = ' ' - fil_loc = ' ' - lwrite_loc = .false. - lread_loc = .false. - iofspin = 1 - ikind = 0 - bdl1 = 0.d0 - bdl2 = 0.d0 - bdr1 = 0.d0 - bdr2 = 0.d0 - nz1 = 11 - dnslab = 0 - energy0 = 0.d0 - denergy = 0.d0 - ecut2d = 0.d0 - ewind = 0.d0 - llocal = .false. - llapack = .false. - epsproj = 1.d-3 - delgep = 0.d0 - cutplot = 2.d0 + outdir = './' + prefix = ' ' + band_file = ' ' + tran_file = ' ' + fil_loc = ' ' + lwrite_loc = .FALSE. + lread_loc = .FALSE. + iofspin = 1 + ikind = 0 + bdl1 = 0.d0 + bdl2 = 0.d0 + bdr1 = 0.d0 + bdr2 = 0.d0 + nz1 = 11 + dnslab = 0 + energy0 = 0.d0 + denergy = 0.d0 + ecut2d = 0.d0 + ewind = 0.d0 + llocal = .FALSE. + llapack = .FALSE. + epsproj = 1.d-3 + delgep = 0.d0 + cutplot = 2.d0 - - -#ifdef __PARA - if (me == 1) then -#endif - ! - ! ... Input from file ? - ! - nargs = iargc() - ! - DO iiarg = 1, ( nargs - 1 ) + IF ( ionode ) THEN ! - CALL getarg( iiarg, input_file ) - IF ( TRIM( input_file ) == '-input' .OR. & - TRIM( input_file ) == '-inp' .OR. & - TRIM( input_file ) == '-in' ) THEN - ! - CALL getarg( ( iiarg + 1 ) , input_file ) - OPEN ( UNIT = 5, FILE = input_file, FORM = 'FORMATTED', & - STATUS = 'OLD', IOSTAT = ierr ) - CALL errore( 'iosys', 'input file ' // TRIM( input_file ) // & - & ' not found' , ierr ) - ! - END IF + ! ... Input from file ? ! - END DO + nargs = iargc() + ! + DO iiarg = 1, ( nargs - 1 ) + ! + CALL getarg( iiarg, input_file ) + IF ( TRIM( input_file ) == '-input' .OR. & + TRIM( input_file ) == '-inp' .OR. & + TRIM( input_file ) == '-in' ) THEN + ! + CALL getarg( ( iiarg + 1 ) , input_file ) + OPEN ( UNIT = 5, FILE = input_file, FORM = 'FORMATTED', & + STATUS = 'OLD', IOSTAT = ierr ) + CALL errore( 'iosys', 'input file ' // TRIM( input_file ) // & + & ' not found' , ierr ) + ! + END IF + ! + END DO + ! + ! reading the namelist inputpp + ! + READ (5, inputcond, err=200, iostat=ios ) +200 CALL errore ('do_cond','reading inputcond namelist',ABS(ios)) + tmp_dir=TRIM(outdir) + ! + ! Reading 2D k-point + READ(5, *, err=300, iostat=ios ) nkpts + ALLOCATE( xyk(2,nkpts) ) + ALLOCATE( wkpt(nkpts) ) + wtot = 0.d0 + DO ik = 1, nkpts + READ(5, *, err=300, iostat=ios) xyk(1,ik), xyk(2,ik), wkpt(ik) + wtot = wtot + wkpt(ik) + ENDDO + DO ik = 1, nkpts + wkpt(ik) = wkpt(ik)/wtot + ENDDO +300 CALL errore ('do_cond','2D k-point',ABS(ios)) - -! -! reading the namelist inputpp -! - read (5, inputcond, err=200, iostat=ios ) -200 call errore ('do_cond','reading inputcond namelist',abs(ios)) - tmp_dir=trim(outdir) -! -! Reading 2D k-point - read(5, *, err=300, iostat=ios ) nkpts - allocate( xyk(2,nkpts) ) - allocate( wkpt(nkpts) ) - wtot = 0.d0 - do ik = 1, nkpts - read(5, *, err=300, iostat=ios) xyk(1,ik), xyk(2,ik), wkpt(ik) - wtot = wtot + wkpt(ik) - enddo - do ik = 1, nkpts - wkpt(ik) = wkpt(ik)/wtot - enddo -300 call errore ('do_cond','2D k-point',abs(ios)) - -! -! To form the array of energies for calculation -! - read(5, *, err=400, iostat=ios ) nenergy - allocate( earr(nenergy) ) - allocate( tran_tot(nenergy) ) - if(abs(denergy).le.1.d-8) then -! the list of energies is read - do ien = 1, nenergy - read(5, *, err=400, iostat=ios) earr(ien) - enddo - else -! the array of energies is automatically formed - do ien = 1, nenergy - earr(ien) = energy0 + (ien-1)*denergy - tran_tot(ien) = 0.d0 - enddo - endif -400 call errore ('do_cond','reading energy list',abs(ios)) -#ifdef __PARA - end if + ! + ! To form the array of energies for calculation + ! + READ(5, *, err=400, iostat=ios ) nenergy + ALLOCATE( earr(nenergy) ) + ALLOCATE( tran_tot(nenergy) ) + IF(ABS(denergy).LE.1.d-8) THEN + ! the list of energies is read + DO ien = 1, nenergy + READ(5, *, err=400, iostat=ios) earr(ien) + ENDDO + ELSE + ! the array of energies is automatically formed + DO ien = 1, nenergy + earr(ien) = energy0 + (ien-1)*denergy + tran_tot(ien) = 0.d0 + ENDDO + ENDIF +400 CALL errore ('do_cond','reading energy list',ABS(ios)) + ! + END IF ! ! ... Broadcast variables @@ -170,88 +165,87 @@ implicit none CALL mp_bcast( llapack, ionode_id ) CALL mp_bcast( nkpts, ionode_id ) CALL mp_bcast( nenergy, ionode_id ) - if (me .ne. 1) then - allocate( xyk(2,nkpts) ) - allocate( wkpt(nkpts) ) - allocate( earr(nenergy) ) - allocate( tran_tot(nenergy) ) - endif + + IF ( .NOT. ionode ) THEN + ALLOCATE( xyk(2,nkpts) ) + ALLOCATE( wkpt(nkpts) ) + ALLOCATE( earr(nenergy) ) + ALLOCATE( tran_tot(nenergy) ) + ENDIF CALL mp_bcast( xyk, ionode_id ) CALL mp_bcast( wkpt, ionode_id ) CALL mp_bcast( earr, ionode_id ) CALL mp_bcast( tran_tot, ionode_id ) -#endif - ! ! Now allocate space for pwscf variables, read and check them. ! - call read_file - call openfil - call struc_fact (nat,tau,ntyp,ityp,ngm,g,bg, & + CALL read_file + CALL openfil + CALL struc_fact (nat,tau,ntyp,ityp,ngm,g,bg, & nr1,nr2,nr3,strf,eigts1,eigts2,eigts3) - call init_us_1 - call newd + CALL init_us_1 + CALL newd ! ! Allocation for pwcond variables ! - call allocate_cond - call init_cond - call stop_clock('init') + CALL allocate_cond + CALL init_cond + CALL stop_clock('init') - if (llocal) & - call local_set(nocrosl,noinsl,noinss,nocrosr,noinsr,norb,norbs) - call poten + IF (llocal) & + CALL local_set(nocrosl,noinsl,noinss,nocrosr,noinsr,norb,norbs) + CALL poten - do ik=1, nkpts + DO ik=1, nkpts - call init_gper(ik) + CALL init_gper(ik) ! ! The main loop ! eryd = earr(1)/rytoev + ef - call local - do ien=1, nenergy + CALL local + DO ien=1, nenergy eryd = earr(ien)/rytoev + ef - call form_zk(n2d, nrzp, zkr, zk, eryd, tpiba) + CALL form_zk(n2d, nrzp, zkr, zk, eryd, tpiba) orbin=1 orbfin=orbin-1+2*nocrosl+noinsl - call scatter_forw(bdl1, bdl2, orbin, orbfin) + CALL scatter_forw(bdl1, bdl2, orbin, orbfin) orbin=1 orbfin=orbin-1+2*nocrosl+noinsl - call compbs(0, bdl1, bdl2, nocrosl, & + CALL compbs(0, bdl1, bdl2, nocrosl, & orbfin-orbin+1, orbin, nchanl, kvall, & kfunl, kfundl, kintl, kcoefl) - if (ikind.eq.2) then + IF (ikind.EQ.2) THEN orbin=2*nocrosl+noinsl+noinss+1 orbfin=orbin-1+2*nocrosr+noinsr - call scatter_forw(bdr1, bdr2, orbin, orbfin) + CALL scatter_forw(bdr1, bdr2, orbin, orbfin) orbin=2*nocrosl+noinsl+noinss+1 orbfin=orbin-1+2*nocrosr+noinsr - call compbs(1, bdr1, bdr2, nocrosr, & + CALL compbs(1, bdr1, bdr2, nocrosr, & orbfin-orbin+1, orbin, nchanr, kvalr, & kfunr, kfundr, kintr, kcoefr) - endif - call summary_band(ik,ien) - if (ikind.ne.0) then + ENDIF + CALL summary_band(ik,ien) + IF (ikind.NE.0) THEN orbin=nocrosl+noinsl+1 orbfin=orbin-1+nocrosl+noinss+nocrosr - call scatter_forw(bdl2, bdr1, orbin, orbfin) - call transmit(ik,ien) - endif - enddo - call free_mem + CALL scatter_forw(bdl2, bdr1, orbin, orbfin) + CALL transmit(ik,ien) + ENDIF + ENDDO + CALL free_mem - enddo + ENDDO - if(ikind.gt.0.and.tran_file.ne.' ') & - call summary_tran(tran_file,nenergy,earr,tran_tot) + IF(ikind.GT.0.AND.tran_file.NE.' ') & + CALL summary_tran(tran_file,nenergy,earr,tran_tot) - call print_clock_pwcond() - call stop_clock('PWCOND') - return -end subroutine do_cond + CALL print_clock_pwcond() + CALL stop_clock('PWCOND') + RETURN +END SUBROUTINE do_cond diff --git a/PWCOND/local.f90 b/PWCOND/local.f90 index 98d2f843b..62175b6ab 100644 --- a/PWCOND/local.f90 +++ b/PWCOND/local.f90 @@ -1,4 +1,3 @@ - ! ! Copyright (C) 2003 A. Smogunov ! This file is distributed under the terms of the @@ -8,100 +7,97 @@ ! ! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC). ! +#include "f_defs.h" ! -subroutine local +SUBROUTINE local ! ! This subroutine computes 2D eigenfunctions and eigenvalues for ! the local potential in each slab and performs 2D reduction of ! the plane wave basis set. ! -#include "f_defs.h" - USE io_global, ONLY : stdout + + USE io_global, ONLY : stdout, ionode USE pwcom USE noncollin_module, ONLY : npol USE io_files USE cond -#ifdef __PARA - USE mp_global, ONLY: nproc - use para + USE mp_global, ONLY : nproc, me_pool, root_pool USE parallel_include -#endif - implicit none + IMPLICIT NONE - integer :: i, il, j, jl, ixy, ig, jg, ipol, igper, k, & + INTEGER :: i, il, j, jl, ixy, ig, jg, ipol, igper, k, & ios, index, number, nprob, nteam, nteamnow, & status, info, kin, kfin, is, js - integer, allocatable :: fftxy(:,:) - real(kind=DP), parameter :: eps=1.d-6 - real(kind=DP), allocatable :: el(:), gp(:) - complex(kind=DP), allocatable :: amat(:,:), amat1(:,:), ymat(:,:), & + INTEGER, ALLOCATABLE :: fftxy(:,:) + REAL(kind=DP), PARAMETER :: eps=1.d-6 + REAL(kind=DP), ALLOCATABLE :: el(:), gp(:) + COMPLEX(kind=DP), ALLOCATABLE :: amat(:,:), amat1(:,:), ymat(:,:), & psibase(:,:), psiprob(:,:) - complex(kind=DP),parameter :: one=(1.d0,0.d0), zero=(0.d0,0.d0) - complex(kind=DP) :: aij, xfact, ZDOTC - logical :: exst + COMPLEX(kind=DP),PARAMETER :: one=(1.d0,0.d0), zero=(0.d0,0.d0) + COMPLEX(kind=DP) :: aij, xfact, ZDOTC + LOGICAL :: exst ! ! To divide the slabs between CPU ! - call start_clock('local') - call slabcpu(nrz, nrzp, nkofz, bdl1, bdl2, bdr1, bdr2, z) + CALL start_clock('local') + CALL slabcpu(nrz, nrzp, nkofz, bdl1, bdl2, bdr1, bdr2, z) ! ! If all the information is already contained in the file it reads it. ! - if (lread_loc) then - call seqopn(4,fil_loc,'unformatted',exst) - if(.not.exst) call errore ('local','fil_loc not found',1) - read(4) n2d + IF (lread_loc) THEN + CALL seqopn(4,fil_loc,'unformatted',exst) + IF(.NOT.exst) CALL errore ('local','fil_loc not found',1) + READ(4) n2d ! Allocate variables depending on n2d - call allocate_cond_2 - read(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d) + CALL allocate_cond_2 + READ(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d) ! WRITE( stdout,*) 'ngper, n2d = ', ngper, n2d - read(4) (((psiper(ig,il,k),ig=1,n2d),il=1,n2d), & + READ(4) (((psiper(ig,il,k),ig=1,n2d),il=1,n2d), & k=1,nrzp) - read(4) ((zkr(il,k),il=1,n2d),k=1,nrzp) + READ(4) ((zkr(il,k),il=1,n2d),k=1,nrzp) - close(unit=4) - return - endif + CLOSE(unit=4) + RETURN + ENDIF - allocate( gp( 2 ) ) - allocate( el( ngper * npol ) ) - allocate( amat( ngper * npol, ngper * npol ) ) - allocate( psibase( ngper * npol, ngper * npol ) ) - allocate( psiprob( ngper * npol, ngper * npol ) ) - allocate( fftxy(-(nrx-1)/2:nrx/2,-(nry-1)/2:nry/2) ) + ALLOCATE( gp( 2 ) ) + ALLOCATE( el( ngper * npol ) ) + ALLOCATE( amat( ngper * npol, ngper * npol ) ) + ALLOCATE( psibase( ngper * npol, ngper * npol ) ) + ALLOCATE( psiprob( ngper * npol, ngper * npol ) ) + ALLOCATE( fftxy(-(nrx-1)/2:nrx/2,-(nry-1)/2:nry/2) ) ! ! To form fftxy correspondence ! - do i=1, nrx + DO i=1, nrx il=i-1 - if (il.gt.nrx/2) il=il-nrx - do j=1, nry + IF (il.GT.nrx/2) il=il-nrx + DO j=1, nry jl=j-1 - if (jl.gt.nry/2) jl=jl-nry + IF (jl.GT.nry/2) jl=jl-nry fftxy(il,jl)=i+(j-1)*nrx - enddo - enddo + ENDDO + ENDDO ! ! to find kin and kfin ! - do k=1, nrz - if (z(k).le.bdl1+eps) kin=k - if (z(k).le.bdr2-eps) kfin=k - enddo + DO k=1, nrz + IF (z(k).LE.bdl1+eps) kin=k + IF (z(k).LE.bdr2-eps) kfin=k + ENDDO ! ! Starting k and number of CPU ! nteam=1 -#ifdef __PARA - kin=kin+me-1 - nteam=nproc -#endif + + kin = kin + me_pool + nteam = nproc ! ! set and solve the eigenvalue equation for each slab @@ -110,31 +106,31 @@ subroutine local nprob=0 psibase=(0.d0,0.d0) - do while(kin.le.kfin) + DO WHILE(kin.LE.kfin) amat=(0.d0,0.d0) - do ig=1, ngper - do jg=ig, ngper - do ipol=1, 2 + DO ig=1, ngper + DO jg=ig, ngper + DO ipol=1, 2 gp(ipol)=gper(ipol,ig)-gper(ipol,jg) - enddo + ENDDO index=number(gp, at, fftxy, nrx, nry) - if (index.gt.0) then - do is=1,npol - do js=1,npol + IF (index.GT.0) THEN + DO is=1,npol + DO js=1,npol amat(ig+(is-1)*ngper,jg+(js-1)*ngper)=vppot(kin,index,is,js) amat(jg+(js-1)*ngper,ig+(is-1)*ngper)= & DCONJG(amat(ig+(is-1)*ngper,jg+(js-1)*ngper)) - enddo - enddo - endif - enddo - do is=1,npol + ENDDO + ENDDO + ENDIF + ENDDO + DO is=1,npol amat(ig+(is-1)*ngper,ig+(is-1)*ngper)= & amat(ig+(is-1)*ngper,ig+(is-1)*ngper)+ & (gper(1,ig)**2 + gper(2,ig)**2)*tpiba2 - enddo - enddo - call hev_ab(ngper*npol, amat, ngper*npol, el, psiprob, & + ENDDO + ENDDO + CALL hev_ab(ngper*npol, amat, ngper*npol, el, psiprob, & -1.d1, eryd+ewind, nprob) ! do is=1,ngper*npol @@ -147,173 +143,173 @@ subroutine local ! stop #ifdef __PARA - if ( me.ne.1 ) then - call mpi_send(nprob,1,MPI_INTEGER,0,17, & + IF ( me_pool == root_pool ) THEN + CALL mpi_send(nprob,1,MPI_INTEGER,0,17, & MPI_COMM_WORLD,info ) - call errore ('n2d reduction','info<>0 in send',info) - call mpi_send(psiprob,2*ngper*npol*ngper*npol,MPI_REAL8,0,18, & + CALL errore ('n2d reduction','info<>0 in send',info) + CALL mpi_send(psiprob,2*ngper*npol*ngper*npol,MPI_REAL8,0,18, & MPI_COMM_WORLD,info ) - call errore ('n2d reduction','info<>0 in send',info) - else - call gramsh(ngper*npol,nprob,1,nprob, & + CALL errore ('n2d reduction','info<>0 in send',info) + ELSE + CALL gramsh(ngper*npol,nprob,1,nprob, & psibase,psiprob,n2d,epsproj) nteamnow=kfin-kin+1 - if(nteamnow.gt.nteam) nteamnow=nteam + IF(nteamnow.GT.nteam) nteamnow=nteam - do ig=1, nteamnow-1 - call mpi_recv(nprob,1,MPI_INTEGER, & + DO ig=1, nteamnow-1 + CALL mpi_recv(nprob,1,MPI_INTEGER, & ig,17,MPI_COMM_WORLD,status,info ) - call errore ('n2d reduction','info<>0 in recv',info) - call mpi_recv(psiprob,2*ngper*npol*ngper*npol,MPI_REAL8, & + CALL errore ('n2d reduction','info<>0 in recv',info) + CALL mpi_recv(psiprob,2*ngper*npol*ngper*npol,MPI_REAL8, & ig,18,MPI_COMM_WORLD,status,info ) - call errore ('n2d reduction','info<>0 in recv',info) - call gramsh(ngper*npol,nprob,1,nprob, & + CALL errore ('n2d reduction','info<>0 in recv',info) + CALL gramsh(ngper*npol,nprob,1,nprob, & psibase,psiprob,n2d,epsproj) - enddo - endif + ENDDO + ENDIF #else - call gramsh(ngper*npol,nprob,1,nprob,psibase,psiprob,n2d,epsproj) + CALL gramsh(ngper*npol,nprob,1,nprob,psibase,psiprob,n2d,epsproj) #endif kin=kin+nteam - enddo + ENDDO #ifdef __PARA - call mpi_barrier( MPI_COMM_WORLD, info ) - call mpi_bcast(n2d,1,MPI_INTEGER,0,MPI_COMM_WORLD,info) - call errore ('reduction','mpi_bcast 1',info) - call mpi_bcast(psibase,2*ngper*npol*ngper*npol,MPI_REAL8,0, & + CALL mpi_barrier( ) + CALL mpi_bcast(n2d,1,MPI_INTEGER,0,MPI_COMM_WORLD,info) + CALL errore ('reduction','mpi_bcast 1',info) + CALL mpi_bcast(psibase,2*ngper*npol*ngper*npol,MPI_REAL8,0, & MPI_COMM_WORLD,info) - call errore ('reduction','mpi_bcast 1',info) + CALL errore ('reduction','mpi_bcast 1',info) #endif ! ! Allocate variables depending on n2d ! - call allocate_cond_2 - if (npol.eq.2) then + CALL allocate_cond_2 + IF (npol.EQ.2) THEN WRITE( stdout,*) 'ngper, ngper*npol, n2d = ', ngper, ngper*npol, n2d - else + ELSE WRITE( stdout,*) 'ngper, n2d = ', ngper, n2d - endif + ENDIF ! ! Construct components of basis vector set on G_per ! - call DCOPY(2*ngper*npol*n2d,psibase,1,newbg,1) + CALL DCOPY(2*ngper*npol*n2d,psibase,1,newbg,1) ! ! set and solve the eigenvalue equation for each slab ! - allocate( amat1( n2d, n2d ) ) - allocate( ymat( ngper*npol, n2d ) ) + ALLOCATE( amat1( n2d, n2d ) ) + ALLOCATE( ymat( ngper*npol, n2d ) ) ! for reduced basis set H'_{ab}=e*^i_aH_{ij}e^j_b - do k=1, nrz - if(nkofz(k).ne.0) then + DO k=1, nrz + IF(nkofz(k).NE.0) THEN ymat=(0.d0,0.d0) ! ! First compute y_{ib}=H_{ij}e_{jb} ! - do ig=1, ngper - do jg=1, ngper - do ipol=1, 2 + DO ig=1, ngper + DO jg=1, ngper + DO ipol=1, 2 gp(ipol) = gper(ipol,ig) - gper(ipol,jg) - enddo + ENDDO index=number(gp, at, fftxy, nrx, nry) - do is=1,npol - do js=1,npol - if (index.gt.0) then + DO is=1,npol + DO js=1,npol + IF (index.GT.0) THEN aij=vppot(k,index,is,js) - else + ELSE aij=(0.d0,0.d0) - endif - if ((ig.eq.jg).and.(is.eq.js)) & + ENDIF + IF ((ig.EQ.jg).AND.(is.EQ.js)) & aij=aij+(gper(1,ig)**2+ & gper(2,ig)**2)*tpiba2 amat(ig+(is-1)*ngper,jg+(js-1)*ngper)= aij - enddo - enddo - enddo - enddo - call ZGEMM('n','n',ngper*npol,n2d,ngper*npol,one,amat,ngper*npol, & + ENDDO + ENDDO + ENDDO + ENDDO + CALL ZGEMM('n','n',ngper*npol,n2d,ngper*npol,one,amat,ngper*npol, & newbg,ngper*npol,zero,ymat,ngper*npol) ! ! and construct H'_{ab}= ! - do il=1, n2d - do jl=il, n2d + DO il=1, n2d + DO jl=il, n2d amat1(il,jl)=ZDOTC(ngper*npol,newbg(1,il),1,ymat(1,jl),1) - amat1(jl,il)=conjg(amat1(il,jl)) - enddo - enddo + amat1(jl,il)=CONJG(amat1(il,jl)) + ENDDO + ENDDO ! ! Solving the eigenvalue problem and construction zk ! info=-1 - call hev_ab(n2d, amat1, n2d, zkr(1,nkofz(k)), & + CALL hev_ab(n2d, amat1, n2d, zkr(1,nkofz(k)), & psiper(1,1,nkofz(k)), 0.d0, 0.d0, info) - endif - enddo + ENDIF + ENDDO #ifdef __PARA - call mpi_barrier( MPI_COMM_WORLD, info ) + CALL mpi_barrier() #endif ! ! saving the 2D data on the file if lwrite_loc=.t. ! - if (lwrite_loc) then - if(fil_loc.eq.' ') call errore ('local','fil_loc no name',1) - call seqopn(4,fil_loc,'unformatted',exst) - write(4) n2d - write(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d) + IF (lwrite_loc) THEN + IF(fil_loc.EQ.' ') CALL errore ('local','fil_loc no name',1) + CALL seqopn(4,fil_loc,'unformatted',exst) + WRITE(4) n2d + WRITE(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d) - write(4) (((psiper(ig,il,k),ig=1,n2d),il=1,n2d), & + WRITE(4) (((psiper(ig,il,k),ig=1,n2d),il=1,n2d), & k=1,nrzp) - write(4) ((zkr(il,k),il=1,n2d),k=1,nrzp) - close(unit=4) - endif + WRITE(4) ((zkr(il,k),il=1,n2d),k=1,nrzp) + CLOSE(unit=4) + ENDIF - deallocate(amat) - deallocate(amat1) - deallocate(ymat) - deallocate(gp) - deallocate(psibase) - deallocate(psiprob) - deallocate(el) - deallocate(fftxy) + DEALLOCATE(amat) + DEALLOCATE(amat1) + DEALLOCATE(ymat) + DEALLOCATE(gp) + DEALLOCATE(psibase) + DEALLOCATE(psiprob) + DEALLOCATE(el) + DEALLOCATE(fftxy) - call stop_clock('local') - return -end subroutine local + CALL stop_clock('local') + RETURN +END SUBROUTINE local !----------------------------------- -function number(gp, at, fftxy, nrx, nry) +FUNCTION number(gp, at, fftxy, nrx, nry) ! ! This function receives as input the coordinates of 2D g vector ! and write on output its fft position. ! - implicit none - integer :: nrx, nry, fftxy(-(nrx-1)/2:nrx/2, -(nry-1)/2:nry/2), & + IMPLICIT NONE + INTEGER :: nrx, nry, fftxy(-(nrx-1)/2:nrx/2, -(nry-1)/2:nry/2), & number, n1, n2 - real(kind=kind(0.d0)) :: gp(2), at(3,3), x1, x2 - real(kind=kind(0.d0)), parameter :: eps=1.d-4 + REAL(kind=KIND(0.d0)) :: gp(2), at(3,3), x1, x2 + REAL(kind=KIND(0.d0)), PARAMETER :: eps=1.d-4 x1=gp(1)*at(1,1)+gp(2)*at(2,1) x2=gp(1)*at(1,2)+gp(2)*at(2,2) - n1=nint(x1) - n2=nint(x2) - if (n1.le.nrx/2.and.n1.ge.-(nrx-1)/2.and. & - n2.le.nry/2.and.n2.ge.-(nry-1)/2) then + n1=NINT(x1) + n2=NINT(x2) + IF (n1.LE.nrx/2.AND.n1.GE.-(nrx-1)/2.AND. & + n2.LE.nry/2.AND.n2.GE.-(nry-1)/2) THEN number=fftxy(n1,n2) - else + ELSE ! ! The g vector is not inside the 2D mesh ! number=-1 - endif + ENDIF - return -end function number + RETURN +END FUNCTION number diff --git a/PWCOND/poten.f90 b/PWCOND/poten.f90 index 5bda405c9..632a0c86c 100644 --- a/PWCOND/poten.f90 +++ b/PWCOND/poten.f90 @@ -8,131 +8,131 @@ ! ! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC). ! +#include "f_defs.h" ! -subroutine poten +SUBROUTINE poten ! ! This subroutine computes the 2D Fourier components of the ! local potential in each slab. ! -#include "f_defs.h" - use pwcom - use noncollin_module, ONLY : noncolin, npol - use cond -#ifdef __PARA - use para -#endif - implicit none + USE pwcom + USE noncollin_module, ONLY : noncolin, npol + USE cond + USE mp_global, ONLY : me_pool + USE pfft, ONLY : npp - integer :: & + IMPLICIT NONE + + INTEGER :: & i, j, ij, ijx, k, n, p, il, ik, kstart, klast, & ix, jx, kx, ir, ir1, ixy, info - integer :: iis, jjs, is(4), js(4), ispin, nspin_eff - integer :: ionode_id - integer, allocatable :: ipiv(:) + INTEGER :: iis, jjs, is(4), js(4), ispin, nspin_eff + INTEGER :: ionode_id + INTEGER, ALLOCATABLE :: ipiv(:) - real(kind=DP), parameter :: eps = 1.d-8 - real(kind=DP) :: arg, bet - real(kind=DP), allocatable :: gz(:), auxr(:) + REAL(kind=DP), PARAMETER :: eps = 1.d-8 + REAL(kind=DP) :: arg, bet + REAL(kind=DP), ALLOCATABLE :: gz(:), auxr(:) - complex(kind=DP), parameter :: cim = (0.d0,1.d0) - complex(kind=DP) :: caux - complex(kind=DP), allocatable :: aux(:), amat(:,:), amat0(:,:) - complex(kind=DP), allocatable :: vppot0(:,:,:,:) + COMPLEX(kind=DP), PARAMETER :: cim = (0.d0,1.d0) + COMPLEX(kind=DP) :: caux + COMPLEX(kind=DP), ALLOCATABLE :: aux(:), amat(:,:), amat0(:,:) + COMPLEX(kind=DP), ALLOCATABLE :: vppot0(:,:,:,:) - logical :: lg + LOGICAL :: lg - call start_clock('poten') - allocate( ipiv( nrz ) ) - allocate( gz( nrz ) ) - allocate( aux( nrx1*nrx2*nrx3 ) ) - allocate( auxr( nrxx ) ) - allocate( amat( nrz, nrz ) ) - allocate( amat0( nrz, nrz ) ) + CALL start_clock('poten') + ALLOCATE( ipiv( nrz ) ) + ALLOCATE( gz( nrz ) ) + ALLOCATE( aux( nrx1*nrx2*nrx3 ) ) + ALLOCATE( auxr( nrxx ) ) + ALLOCATE( amat( nrz, nrz ) ) + ALLOCATE( amat0( nrz, nrz ) ) ! ! Compute the Gz vectors in the z direction ! - do k = 1, nrz + DO k = 1, nrz il = k-1 - if (il.gt.nrz/2) il = il-nrz + IF (il.GT.nrz/2) il = il-nrz gz(k) = il*bg(3,3) - enddo + ENDDO ! ! set up the matrix for the linear system ! -do n=1,nrz - do p=1,nrz +DO n=1,nrz + DO p=1,nrz arg=gz(n)*z(p)*tpi bet=gz(n)*(z(p+1)-z(p))*tpi - if (abs(gz(n)).gt.eps) then - caux=cim*(CMPLX(cos(bet),-sin(bet))-(1.d0,0.d0)) & + IF (ABS(gz(n)).GT.eps) THEN + caux=cim*(CMPLX(COS(bet),-SIN(bet))-(1.d0,0.d0)) & /zl/gz(n)/tpi - else + ELSE caux=(z(p+1)-z(p))/zl - endif - amat0(n,p)=CMPLX(cos(arg),-sin(arg))*caux - enddo -enddo -if (noncolin) then + ENDIF + amat0(n,p)=CMPLX(COS(arg),-SIN(arg))*caux + ENDDO +ENDDO +IF (noncolin) THEN nspin_eff=4 ij=0 - do iis=1,2 - do jjs=1,2 + DO iis=1,2 + DO jjs=1,2 ij=ij+1 is(ij)=iis js(ij)=jjs - enddo - enddo -else + ENDDO + ENDDO +ELSE nspin_eff=1 is(1)=1 js(1)=1 -endif +ENDIF ! ! To form local potential on the real space mesh ! vppot = 0.d0 -do ispin=1,nspin_eff - if (noncolin) then - if (ispin==1) then +DO ispin=1,nspin_eff + IF (noncolin) THEN + IF (ispin==1) THEN auxr(:) = vltot(:)+vr(:,1) - else + ELSE auxr(:) = vr(:,ispin) - endif - else + ENDIF + ELSE auxr(:) = vltot(:) + vr(:,iofspin) - endif + ENDIF ! ! To collect the potential from different CPUs ! aux(:) = (0.d0,0.d0) #ifdef __PARA kstart = 1 - do i=1, me-1 + DO i=1, me_pool kstart = kstart+npp(i) - enddo - klast = kstart+npp(me)-1 + ENDDO + klast = kstart+npp(me_pool+1)-1 #endif - do i = 1, nrx1*nrx2 - do k = 1, nr3 - lg = .true. + DO i = 1, nrx1*nrx2 + DO k = 1, nr3 + lg = .TRUE. ir = i+(k-1)*nrx2*nrx1 ir1 = ir #ifdef __PARA - if(k.ge.kstart.and.k.le.klast) then - lg = .true. + IF(k.GE.kstart.AND.k.LE.klast) THEN + lg = .TRUE. ir1 = i+(k-kstart)*nrx2*nrx1 - else - lg = .false. - endif + ELSE + lg = .FALSE. + ENDIF #endif - if (lg) then + IF (lg) THEN aux(ir) = auxr(ir1) - endif - enddo - enddo + ENDIF + ENDDO + ENDDO #ifdef __PARA - call reduce (2*nrx1*nrx2*nrx3,aux) + CALL reduce (2*nrx1*nrx2*nrx3,aux) #endif @@ -143,74 +143,74 @@ do ispin=1,nspin_eff ! ! This FFT is needed to make a non-parallel FFT in the parallel case ! - call cft3sp(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1) + CALL cft3sp(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1) #else - call cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1) + CALL cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1) #endif - do i = 1, nrx - if(i.gt.nrx/2+1) then + DO i = 1, nrx + IF(i.GT.nrx/2+1) THEN ix = nr1-(nrx-i) - else + ELSE ix = i - endif - do j = 1, nry - if(j.gt.nry/2+1) then + ENDIF + DO j = 1, nry + IF(j.GT.nry/2+1) THEN jx = nr2-(nry-j) - else + ELSE jx = j - endif + ENDIF ij = i+(j-1)*nrx ijx = ix+(jx-1)*nrx1 - do k = 1, nrz + DO k = 1, nrz il = k-1 - if (il.gt.nrz/2) il = il-nrz - if(il.le.nr3/2.and.il.ge.-(nr3-1)/2) then + IF (il.GT.nrz/2) il = il-nrz + IF(il.LE.nr3/2.AND.il.GE.-(nr3-1)/2) THEN - if(k.gt.nrz/2+1) then + IF(k.GT.nrz/2+1) THEN kx = nr3-(nrz-k) - else + ELSE kx = k - endif + ENDIF vppot(k, ij, is(ispin), js(ispin)) = aux(ijx+(kx-1)*nrx1*nrx2) - endif - enddo - enddo - enddo + ENDIF + ENDDO + ENDDO + ENDDO ! ! solve the linear system ! amat=amat0 - call ZGESV(nrz, nrx*nry, amat, nrz, ipiv, vppot(1,1,is(ispin),js(ispin)),& + CALL ZGESV(nrz, nrx*nry, amat, nrz, ipiv, vppot(1,1,is(ispin),js(ispin)),& nrz, info) - call errore ('poten','info different from zero',abs(info)) -enddo + CALL errore ('poten','info different from zero',ABS(info)) +ENDDO -if (noncolin) then - allocate( vppot0(nrz, nrx * nry, npol, npol) ) +IF (noncolin) THEN + ALLOCATE( vppot0(nrz, nrx * nry, npol, npol) ) vppot0=vppot vppot(:,:,1,1)=vppot0(:,:,1,1)+vppot0(:,:,2,2) vppot(:,:,1,2)=vppot0(:,:,1,2)-(0.d0,1.d0)*vppot0(:,:,2,1) vppot(:,:,2,1)=vppot0(:,:,1,2)+(0.d0,1.d0)*vppot0(:,:,2,1) vppot(:,:,2,2)=vppot0(:,:,1,1)-vppot0(:,:,2,2) - deallocate( vppot0 ) -endif + DEALLOCATE( vppot0 ) +ENDIF ! do p = 1, nrz ! write(6,'(i5,2f12.6)') p, real(vppot(p,105,2,2)), imag(vppot(p,105,2,2)) ! enddo ! stop - deallocate(ipiv) - deallocate(gz) - deallocate(aux) - deallocate(auxr) - deallocate(amat) - deallocate(amat0) + DEALLOCATE(ipiv) + DEALLOCATE(gz) + DEALLOCATE(aux) + DEALLOCATE(auxr) + DEALLOCATE(amat) + DEALLOCATE(amat0) - call stop_clock('poten') + CALL stop_clock('poten') - return -end subroutine poten + RETURN +END SUBROUTINE poten diff --git a/PWCOND/scatter_forw.f90 b/PWCOND/scatter_forw.f90 index c639492aa..6468168a1 100644 --- a/PWCOND/scatter_forw.f90 +++ b/PWCOND/scatter_forw.f90 @@ -9,8 +9,9 @@ ! Optimized Aug. 2004 (ADC) ! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC). ! +#include "f_defs.h" ! -subroutine scatter_forw(zin, zfin, orbin, orbfin) +SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin) ! ! This subroutine computes local Phi_n and partial nonlocal Phi_alpha ! solutions of the Schrodinger equation in the region zin