1c
2c     $Id$
3c
4
5*  ************************************************************
6*  *                                                          *
7*  *             Band by Band Kohn-Sham Minimizer             *
8*  *                                                          *
9*  *                                                          *
10*  *                                                          *
11*  ************************************************************
12
13      subroutine bybminimize(E,deltae,deltac,current_iteration,
14     >                       set_iterations,iterations,failed)
15      implicit none
16      real*8     E(*)
17      real*8     deltae,deltac
18      integer    current_iteration
19      logical    set_iterations
20      integer    iterations
21      logical    failed
22
23
24#include "stdio.fh"
25#include "bafdecls.fh"
26#include "util.fh"
27
28*     **** local variables ****
29      integer MAX_SD_COUNT
30      parameter (MAX_SD_COUNT = 3)
31      integer MASTER,taskid
32      parameter (MASTER=0)
33
34      real*8  deltat_min
35      parameter (deltat_min=1.0d-3)
36
37      integer vall_in(2),vall_out(2),vall_junk(2),rho_in(2)
38      real*8  E0,dE0,deltae_old,Ein,deltae_history(10)
39      real*8  ks_deltae,deltav,dV,deltav_old,diis_error
40      integer nx,ny,nz,stalled_count,sd_count
41
42
43      real*8     tole,tolc
44      real*8     ehartree,eorbit,exc,pxc,eion
45      real*8     Enew,Eold,alpha
46      !common / cgsd_block / Enew,Eold,alpha
47
48      integer it,it_in,i,ispin,bfgscount,icount,sd_it,cg_it
49
50      logical value,precondition,done,stalled,deltav_bad(4),oprint
51      integer n2ft3d,n2ft3d_map
52      !real*8  e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
53      !real*8  e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
54      real*8 e_lj,e_q,e_spring
55      real*8 ehfx,phfx
56
57
58
59*     **** external functions ****
60      logical control_print
61      integer  control_ispin,control_scf_algorithm,control_ks_algorithm
62      integer  control_it_in,control_it_out,psi_ne,control_version
63      real*8   control_tole,control_tolc,control_ks_alpha
64      real*8   rho_error,psi_1energy,psi_error
65      real*8   dng_1ehartree,lattice_omega
66      real*8   psi_1ke
67      real*8   psi_1vl,psi_1v_field,dng_1vl_mm
68      real*8   psi_1vnl
69      real*8   rho_1exc
70      real*8   rho_1pxc
71      real*8   ewald_e,ion_ion_e
72      real*8   psi_1eorbit
73
74      external control_print
75      external control_ispin,control_scf_algorithm,control_ks_algorithm
76      external control_it_in,control_it_out,psi_ne,control_version
77      external control_tole,control_tolc,control_ks_alpha
78      external rho_error,psi_1energy,psi_error
79      external dng_1ehartree,lattice_omega
80      external psi_1ke
81      external psi_1vl,psi_1v_field,dng_1vl_mm
82      external psi_1vnl
83      external rho_1exc
84      external rho_1pxc
85      external ewald_e,ion_ion_e
86      external psi_1eorbit
87
88*     ***** QM/MM external functions ****
89      logical  pspw_qmmm_found
90      real*8   pspw_qmmm_LJ_E,pspw_qmmm_Q_E,pspw_qmmm_spring_E
91      real*8   pspw_qmmm_LJ_Emix,pspw_qmmm_Q_Emix
92      external pspw_qmmm_found
93      external pspw_qmmm_LJ_E,pspw_qmmm_Q_E,pspw_qmmm_spring_E
94      external pspw_qmmm_LJ_Emix,pspw_qmmm_Q_Emix
95
96*     ***** pspw_charge external functions ****
97      logical  pspw_charge_found,control_precondition,pspw_HFX
98      real*8   pspw_charge_Energy_ion,pspw_charge_Energy_charge
99      external pspw_charge_found,control_precondition,pspw_HFX
100      external pspw_charge_Energy_ion,pspw_charge_Energy_charge
101
102      real*8   psi_1_noupdate_energy,psi_eigenvalue
103      external psi_1_noupdate_energy,psi_eigenvalue
104      logical  psp_U_psputerm,meta_found
105      external psp_U_psputerm,meta_found
106      logical  nwpw_meta_gga_on,ion_disp_on
107      external nwpw_meta_gga_on,ion_disp_on
108      real*8   psi_1meta_gga_pxc,ion_disp_energy
109      external psi_1meta_gga_pxc,ion_disp_energy
110
111
112
113      Ein = E(1)
114      call Parallel_taskid(taskid)
115      oprint = (taskid.eq.MASTER).and.control_print(print_medium)
116
117      call D3dB_nx(1,nx)
118      call D3dB_ny(1,ny)
119      call D3dB_nz(1,nz)
120      dV = lattice_omega()/dble(nx*ny*nz)
121      if (set_iterations) then
122        it_in = iterations
123        sd_it = 2
124        cg_it = 1
125      else
126        it_in = control_it_in()*control_it_out()
127        sd_it = 10
128        cg_it = 10
129      end if
130      tole  = control_tole()
131      tolc  = control_tolc()
132      precondition = control_precondition()
133      ispin = control_ispin()
134      deltav_old = 10.0d0
135      deltav     = 0.0d0
136
137      stalled       = .false.
138      deltae_history(1) = 0.0d0
139      deltae_history(2) = 0.0d0
140      deltae_history(3) = 0.0d0
141      deltae_history(4) = 0.0d0
142      stalled_count     = 0
143      sd_count          = 0
144
145      call D3dB_n2ft3d(1,n2ft3d)
146      call D3dB_n2ft3d_map(1,n2ft3d_map)
147
148*     **** allocate rho_in and rho_out ****
149      value = BA_push_get(mt_dbl,2*n2ft3d,
150     >                     'vall_in',vall_in(2),vall_in(1))
151      value = value.and.
152     >        BA_push_get(mt_dbl,2*n2ft3d,
153     >                     'vall_out',vall_out(2),vall_out(1))
154      value = value.and.
155     >        BA_push_get(mt_dbl,2*n2ft3d,
156     >                   'vall_junk',vall_junk(2),vall_junk(1))
157      value = value.and.
158     >        BA_push_get(mt_dbl,2*n2ft3d,
159     >                     'rho_in',rho_in(2),rho_in(1))
160      if (.not. value)
161     >   call errquit('bybminimize:out of stack memory',0,0)
162      call dcopy(2*n2ft3d,0.0d0,0,dbl_mb(vall_in(1)),1)
163      call dcopy(2*n2ft3d,0.0d0,0,dbl_mb(vall_out(1)),1)
164      call dcopy(2*n2ft3d,0.0d0,0,dbl_mb(vall_junk(1)),1)
165      call dcopy(2*n2ft3d,0.0d0,0,dbl_mb(rho_in(1)),1)
166
167*     **** ion-ion energy ****
168      eion = 0.0d0
169      if (control_version().eq.3) eion = ewald_e()
170      if (control_version().eq.4) eion = ion_ion_e()
171
172
173*     **********************
174*     **** bybminimizer ****
175*     **********************
176
177
178*     **** set the initial density ****
179      if (current_iteration.eq.1) then
180         Enew   = psi_1energy()
181         alpha = control_ks_alpha()
182         deltae = -9232323299.0d0
183         ks_deltae = tole
184         call electron_gen_vall()
185         call electron_get_vall(dbl_mb(vall_out(1)))
186         call electron_get_vall(dbl_mb(vall_in(1)))
187
188         call psi_1gen_hml()
189         call psi_diagonalize_hml_assending()
190         call psi_1rotate2()
191         call psi_2to1()
192      else
193         call electron_get_vall(dbl_mb(vall_out(1)))
194         call electron_get_vall(dbl_mb(vall_in(1)))
195         !call psi_get_density(1,dbl_mb(rho_in(1)))
196      end if
197
198*     **** iniitialize SCF Mixing ****
199      call nwpw_scf_mixing_init(control_scf_algorithm(),alpha,
200     >                5,ispin,n2ft3d,dbl_mb(vall_out(1)))
201
202*     **** iniitialize RMM-DIIS ****
203      if (control_ks_algorithm().eq.1) call pspw_rmmdiis_init(5)
204
205
206*     ***** diis loop ****
207      it   = 0
208 2    it   = it + 1
209
210*     **** diaganolize KS matrix ****
211      call psi_KS_update(1,
212     >                   control_ks_algorithm(),
213     >                   precondition,
214     >                   ks_deltae)
215
216      call rho_1to2()
217      Eold = Enew
218      Enew = psi_1energy()
219
220      deltae = Enew-Eold
221
222      call electron_gen_vall()
223      call electron_get_vall(dbl_mb(vall_in(1)))
224
225*     **** compute deltaV ****
226      call dcopy(ispin*n2ft3d_map,
227     >           dbl_mb(vall_in(1)),1,
228     >           dbl_mb(vall_junk(1)),1)
229      call daxpy(ispin*n2ft3d_map,
230     >              (-1.0d0),
231     >              dbl_mb(vall_out(1)),1,
232     >              dbl_mb(vall_junk(1)),1)
233      deltav = ddot(ispin*n2ft3d_map,
234     >                 dbl_mb(vall_junk(1)),1,
235     >                 dbl_mb(vall_junk(1)),1)
236      call D3dB_SumAll(deltav)
237      deltav = deltav*dV
238
239
240
241*     **** update vall using density mixing ****
242c      if ((it.le.0) .or.
243c     >    ((dabs(deltae).lt.1.0d1) .and.
244c     >    (deltav       .lt.1.0d1) .and.
245c     >    (.not.stalled          ))) then
246      if ((it.le.0) .or.
247     >    ((deltae.lt.0.0d0) .and.
248     >    (.not.stalled          ))) then
249
250         call nwpw_scf_mixing(dbl_mb(vall_in(1)),dbl_mb(vall_out(1)),
251     >                        deltae,diis_error)
252
253*     **** bad convergence - try fixed step steepest descent ****
254      else
255
256  30     call sdminimize(sd_it)
257         sd_count = sd_count + 1
258         Eold   = Enew
259         Enew   = psi_1energy()
260
261c         if ((Enew.gt.Eold).or.(dabs(Enew-Eold).gt.1.0d-1)) go to 30
262         if ((Enew.gt.Eold).and.(sd_count.lt.MAX_SD_COUNT)) go to 30
263
264         call dcopy(ispin*n2ft3d,
265     >              dbl_mb(vall_out(1)),1,
266     >              dbl_mb(vall_junk(1)),1)
267
268         call electron_gen_vall()
269         call electron_get_vall(dbl_mb(vall_out(1)))
270         call nwpw_scf_mixing_reset(dbl_mb(vall_out(1)))
271
272
273         call daxpy(ispin*n2ft3d,
274     >              (-1.0d0),
275     >              dbl_mb(vall_out(1)),1,
276     >              dbl_mb(vall_junk(1)),1)
277         deltav = ddot(ispin*n2ft3d,
278     >                 dbl_mb(vall_junk(1)),1,
279     >                 dbl_mb(vall_junk(1)),1)
280         call D3dB_SumAll(deltav)
281         deltav = deltav*dV
282
283         call psi_1gen_hml()
284         call psi_diagonalize_hml_assending()
285         call psi_1rotate2()
286         call psi_2to1()
287
288         stalled       = .false.
289         deltae_history(1) = 0.0d0
290         deltae_history(2) = 0.0d0
291         deltae_history(3) = 0.0d0
292         deltae_history(4) = 0.0d0
293         stalled_count     = 0
294
295      end if
296      call electron_set_vall(dbl_mb(vall_out(1)))
297
298
299*     **** tolerance checks ****
300      deltae = Enew-Eold
301      deltac = rho_error()
302      E(1)   = Enew+eion
303
304      if ((oprint).and.(.not.set_iterations)) then
305        write(luout,1310) it,E(1),deltae,deltac,deltav
306        call util_flush(luout)
307      end if
308 1310 FORMAT(I8,E20.10,3E15.5)
309
310
311      !**** set ks_deltae ****
312      ks_deltae = 0.001d0*dabs(deltae)
313      if (ks_deltae.lt.(0.1d0*tole)) ks_deltae = 0.1d0*tole
314      if (ks_deltae.gt.1.0d-4) ks_deltae = 1.0d-4
315      !ks_deltae = 0.1d0*tole
316
317
318
319      deltav_old = deltav
320
321      deltae_history(1)    = deltae_history(2)
322      deltae_history(2)    = deltae_history(3)
323      deltae_history(3)    = deltae_history(4)
324      deltae_history(4)    = deltae
325
326      if (stalled_count .gt.4) then
327        stalled = (deltae_history(4)
328     >            +deltae_history(3)
329     >            +deltae_history(2)
330     >            +deltae_history(1)).gt.0.0d0
331      else
332         stalled = .false.
333      end if
334      stalled_count = stalled_count + 1
335c      stalled = .false.
336      if (deltae.gt.0.0d0) stalled = .true.
337
338      precondition = precondition.and.(dabs(deltae).gt.1*tole)
339
340      done = ( (    (dabs(deltae).lt.tole)
341     >         .and.(deltae.lt.0.0d0)
342     >         .and.(deltac      .lt.tolc))
343     >       .or. (it.ge.it_in)
344     >       .or. (sd_count.ge.MAX_SD_COUNT))
345
346      if (.not.done) go to 2
347
348
349
350*     **** free memory ****
351      call nwpw_scf_mixing_end()
352      if (control_ks_algorithm().eq.1) call pspw_rmmdiis_end()
353      value =           BA_pop_stack(rho_in(2))
354      value = value.and.BA_pop_stack(vall_junk(2))
355      value = value.and.BA_pop_stack(vall_out(2))
356      value = value.and.BA_pop_stack(vall_in(2))
357      if (.not. value)
358     >  call errquit('bybminimize: popping stack',1,0)
359
360c      call psi_check()
361
362
363*     **** set return energies **** - This is duplicate code
364      !Enew     = psi_1energy()
365      eorbit   = psi_1eorbit()
366      ehartree = dng_1ehartree()
367      exc      = rho_1exc()
368      pxc      = rho_1pxc()
369
370      E(1)  = Enew + eion
371      E(2)  = eorbit
372      E(3)  = ehartree
373      E(4)  = exc
374      E(5)  = eion
375      E(6)  = psi_1ke()
376      E(7)  = psi_1vl()
377      E(8)  = psi_1vnl()
378      E(9)  = 2.0d0*ehartree
379      E(10) = pxc
380      if (pspw_qmmm_found()) then
381         e_lj     = pspw_qmmm_LJ_E()
382         e_q      = pspw_qmmm_Q_E()
383         e_spring = pspw_qmmm_spring_E()
384         E(1)  = E(1) + e_lj + e_q + e_spring
385
386         E(11) = e_lj
387         E(12) = e_q
388         E(13) = e_spring
389
390*        **** Eqm-mm energy ***
391         E(14) = pspw_qmmm_LJ_Emix()
392         E(14) = E(14) + pspw_qmmm_Q_Emix()
393         E(14) = E(14) + dng_1vl_mm()
394
395      end if
396
397*     **** get pspw_charge  energies ****
398      if (pspw_charge_found()) then
399         E(19)  = psi_1v_field()
400         E(20)  = pspw_charge_Energy_ion()
401         E(21)  = pspw_charge_Energy_charge()
402         E(1)   = E(1) + E(20) + E(21)
403      end if
404
405*     **** HFX terms ****
406      if (pspw_HFX()) then
407         call electron_HFX_energies(ehfx,phfx)
408         E(26) = ehfx
409         E(27) = phfx
410      end if
411
412*     **** DFT+U terms ****
413      if (psp_U_psputerm()) then
414         call electron_U_energies(ehfx,phfx)
415         E(29) =  ehfx
416         E(30) =  phfx
417      end if
418
419*     **** Metadynamics potential terms ****
420      if (meta_found()) then
421         call electron_meta_energies(ehfx,phfx)
422         E(31) =  ehfx
423         E(32) =  phfx
424      end if
425
426*     **** Metadynamics GGA Tau potential term ****
427      if (nwpw_meta_gga_on()) then
428         E(10) = E(10) - psi_1meta_gga_pxc()
429      end if
430
431*     **** Dispersion energy ****
432      if (ion_disp_on()) then
433         E(33) = ion_disp_energy()
434         E(1)  = E(1) + E(33)
435      end if
436
437
438      failed = (sd_count.ge.MAX_SD_COUNT)
439
440      return
441      end
442
443
444