1 2! Copyright (C) 2017 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6subroutine gndstulr 7use modmain 8use modulr 9use moddftu 10use modmpi 11use modomp 12use modstore 13implicit none 14! local variables 15logical twrite 16integer ik0,ir,lp,nthd 17integer nmix,nwork,n 18real(8) dv 19! allocatable arrays 20integer(8), allocatable :: lock(:) 21real(8), allocatable :: work(:) 22complex(8), allocatable :: evecu(:,:) 23if (xctype(1).lt.0) then 24 write(*,*) 25 write(*,'("Error(gndstulr): ultra long-range does not work with OEP")') 26 write(*,*) 27 stop 28end if 29if (spincore) then 30 write(*,*) 31 write(*,'("Error(gndstulr): ultra long-range does not work with & 32 &spin-polarised cores")') 33 write(*,*) 34 stop 35end if 36! no k-point reduction 37reducek_=reducek 38reducek=0 39! initialise global variables 40call init0 41call init1 42! write the kappa-points to file 43call writekpa 44! write the k+kappa-points to file 45call writekpts 46! read the regular Kohn-Sham potential from file 47call readstate 48! generate the first- and second-variational eigenvectors and eigenvalues for 49! the k+kappa-point set 50call genvsig 51call gencore 52call readfermi 53call linengy 54call genapwlofr 55call gensocfr 56call genevfsv 57call occupy 58! initialise the ultra long-range variables 59call initulr 60if (task.eq.700) then 61! initialise the long-range Kohn-Sham potential and magnetic field 62 call potuinit 63else 64! read in the potential and density from STATE_ULR.OUT 65 call readstulr 66end if 67! initialise the external Coulomb potential 68call vclqinit 69! write the long-range Hamiltonian diagonal blocks to file 70call genhdbulr 71! size of mixing vector (complex array) 72nmix=2*size(vsbsq) 73! determine the size of the mixer work array 74nwork=-1 75call mixerifc(mixtype,nmix,vsbsq,dv,nwork,vsbsq) 76allocate(work(nwork)) 77! initialise the mixer 78iscl=0 79call mixerifc(mixtype,nmix,vsbsq,dv,nwork,work) 80! initialise the OpenMP locks 81allocate(lock(nqpt)) 82do ir=1,nqpt 83 call omp_init_lock(lock(ir)) 84end do 85! set last self-consistent loop flag 86tlast=.false. 87! begin the self-consistent loop 88if (mp_mpi) then 89! open ULR_INFO.OUT file 90 open(60,file='ULR_INFO.OUT',form='FORMATTED') 91! open RMSDVS.OUT 92 open(65,file='RMSDVS.OUT',form='FORMATTED') 93 call writeinfou(60) 94 write(60,*) 95 write(60,'("+------------------------------+")') 96 write(60,'("| Self-consistent loop started |")') 97 write(60,'("+------------------------------+")') 98end if 99do iscl=1,maxscl 100 if (mp_mpi) then 101 write(60,*) 102 write(60,'("+--------------------+")') 103 write(60,'("| Loop number : ",I4," |")') iscl 104 write(60,'("+--------------------+")') 105 end if 106 if (iscl.ge.maxscl) then 107 if (mp_mpi) then 108 write(60,*) 109 write(60,'("Reached self-consistent loops maximum")') 110 end if 111 write(*,*) 112 write(*,'("Warning(gndstulr): failed to reach self-consistency after ",I4,& 113 &" loops")') iscl 114 tlast=.true. 115 end if 116! reset the OpenMP thread variables 117 call omp_reset 118! apply required local operations to the potential and magnetic field 119 call vblocalu 120! zero the density and magnetisation 121 rhormt(:,:,:)=0.d0 122 rhorir(:,:)=0.d0 123 if (spinpol) then 124 magrmt(:,:,:,:)=0.d0 125 magrir(:,:,:)=0.d0 126 end if 127! loop over original k-points 128 call holdthd(nkpt0/np_mpi,nthd) 129!$OMP PARALLEL DEFAULT(SHARED) & 130!$OMP PRIVATE(evecu) & 131!$OMP NUM_THREADS(nthd) 132 allocate(evecu(nstulr,nstulr)) 133!$OMP DO 134 do ik0=1,nkpt0 135! distribute among MPI processes 136 if (mod(ik0-1,np_mpi).ne.lp_mpi) cycle 137! solve the ultra long-range eigenvalue equation 138 call eveqnulr(ik0,evecu) 139! add to the density, magnetisation and current 140 call rhomaguk(ik0,lock,evecu) 141 end do 142!$OMP END DO 143 deallocate(evecu) 144!$OMP END PARALLEL 145 call freethd(nthd) 146 if (np_mpi.gt.1) then 147! broadcast eigenvalue array to every process 148 do ik0=1,nkpt0 149 lp=mod(ik0-1,np_mpi) 150 call mpi_bcast(evalu(:,ik0),nstulr,mpi_double_precision,lp,mpicom,ierror) 151 end do 152! add densities from each process and redistribute 153 n=npcmtmax*natmtot*nqpt 154 call mpi_allreduce(mpi_in_place,rhormt,n,mpi_double_precision,mpi_sum, & 155 mpicom,ierror) 156 n=ngtc*nqpt 157 call mpi_allreduce(mpi_in_place,rhorir,n,mpi_double_precision,mpi_sum, & 158 mpicom,ierror) 159 if (spinpol) then 160 n=npcmtmax*natmtot*ndmag*nqpt 161 call mpi_allreduce(mpi_in_place,magrmt,n,mpi_double_precision,mpi_sum, & 162 mpicom,ierror) 163 n=ngtc*ndmag*nqpt 164 call mpi_allreduce(mpi_in_place,magrir,n,mpi_double_precision,mpi_sum, & 165 mpicom,ierror) 166 end if 167 end if 168! find the occupation numbers and Fermi energy 169 call occupyulr 170! synchronise MPI processes 171 call mpi_barrier(mpicom,ierror) 172! add the core density 173 call rhocoreu 174! perform partial Fourier transform to Q-space 175 call rhomagq 176! determine the muffin-tin and interstitial charges and moments 177 call chargeu 178 call momentu 179! compute the ultra long-range Kohn-Sham potential 180 call potksu 181! mix the old potential and field with the new 182 call mixerifc(mixtype,nmix,vsbsq,dv,nwork,work) 183! multiply the RMS change in potential by the number of Q-points 184 dv=dv*dble(nfqrz) 185! calculate and add the fixed spin moment effective field (after mixing) 186 call fsmbfield 187 call addbfsmu 188 if (mp_mpi) then 189! write eigenvalues to file 190 call writeevalu 191! output energy components 192 call writeengyu(60) 193! output charges 194 call writechg(60) 195! write muffin-tin charges for each R-vector 196 call writechgrmt 197 if (spinpol) then 198! output moments 199 call writemom(60) 200! write muffin-tin moments for each R-vector 201 call writemomrmt 202 end if 203! output effective fields for fixed spin moment calculations 204 if (fsmtype.ne.0) call writefsm(60) 205! check for WRITE file 206 call checkwrite(twrite) 207! write STATE_ULR.OUT file if required 208 if (twrite.or.((nwrite.ge.1).and.(mod(iscl,nwrite).eq.0))) then 209 call writestulr 210 write(60,*) 211 write(60,'("Wrote STATE_ULR.OUT")') 212 end if 213 end if 214! exit self-consistent loop if required 215 if (tlast) goto 10 216! check for convergence 217 if (iscl.ge.2) then 218 if (mp_mpi) then 219 write(60,*) 220 write(60,'("RMS change in Kohn-Sham potential (target) : ",G18.10," (",& 221 &G18.10,")")') dv,epspot 222 flush(60) 223 write(65,'(G18.10)') dv 224 flush(65) 225 end if 226 if (dv.lt.epspot) then 227 if (mp_mpi) then 228 write(60,*) 229 write(60,'("Convergence targets achieved")') 230 end if 231 tlast=.true. 232 end if 233 end if 234! check for STOP file 235 call checkstop 236 if (tstop) tlast=.true. 237! broadcast tlast from master process to all other processes 238 call mpi_bcast(tlast,1,mpi_logical,0,mpicom,ierror) 239! reset the OpenMP thread variables 240 call omp_reset 241end do 24210 continue 243if (mp_mpi) then 244 write(60,*) 245 write(60,'("+------------------------------+")') 246 write(60,'("| Self-consistent loop stopped |")') 247 write(60,'("+------------------------------+")') 248 if (maxscl.gt.1) then 249 call writestulr 250 write(60,*) 251 write(60,'("Wrote STATE_ULR.OUT")') 252 end if 253! close the ULR_INFO.OUT file 254 close(60) 255! close the RMSDVS.OUT file 256 close(65) 257end if 258! destroy the OpenMP locks 259do ir=1,nqpt 260 call omp_destroy_lock(lock(ir)) 261end do 262deallocate(lock,work) 263! restore original parameters 264reducek=reducek_ 265! synchronise MPI processes 266call mpi_barrier(mpicom,ierror) 267end subroutine 268 269