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