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