1! 2! Copyright (C) 2002-2008 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!---------------------------------------------------------------------------- 10SUBROUTINE cprmain( tau_out, fion_out, etot_out ) 11 !---------------------------------------------------------------------------- 12 ! 13 USE kinds, ONLY : DP 14 USE constants, ONLY : bohr_radius_angs, amu_au, au_gpa 15 USE control_flags, ONLY : iprint, isave, thdyn, tpre, iverbosity, & 16 tfor, remove_rigid_rot, taurdr, llondon,& 17 tprnfor, tsdc, lconstrain, lwf, & 18 ndr, ndw, nomore, tsde, textfor, & 19 tortho, tnosee, tnosep, trane, tranp, & 20 tsdp, tcp, tcap, ampre, amprp, tnoseh, & 21 tolp, ortho_eps, ortho_max, & 22 tfirst, tlast !moved here to make 23 !autopilot work 24 USE core, ONLY : rhoc 25 USE uspp_param, ONLY : nhm, nh, ish 26 USE uspp, ONLY : nkb, vkb, becsum, deeq, okvan, nlcc_any 27 USE energies, ONLY : eht, epseu, exc, etot, eself, enl, & 28 ekin, atot, entropy, egrand, enthal, & 29 ekincm, print_energies 30 USE electrons_base, ONLY : nbspx, nbsp, ispin, f, nspin, nbsp_bgrp 31 USE electrons_base, ONLY : nel, iupdwn, nupdwn, nudx, nelt 32 USE electrons_module, ONLY : distribute_c, collect_c 33 USE efield_module, ONLY : efield, epol, tefield, allocate_efield, & 34 efield_update, ipolp, qmat, gqq, evalue,& 35 berry_energy, pberryel, pberryion, & 36 efield2, epol2, tefield2, & 37 allocate_efield2, efield_update2, & 38 ipolp2, qmat2, gqq2, evalue2, & 39 berry_energy2, pberryel2, pberryion2 40 USE ensemble_dft, ONLY : tens, z0t, gibbsfe 41 USE cg_module, ONLY : tcg, cg_update, c0old 42 USE smallbox_gvec, ONLY : ngb 43 USE gvecw, ONLY : ngw 44 USE gvect, ONLY : gstart, mill, eigts1, eigts2, eigts3 45 USE ions_base, ONLY : na, nat, amass, nax, nsp, rcmax 46 USE ions_base, ONLY : ions_cofmass, ions_kinene, & 47 ions_temp, ions_thermal_stress, & 48 if_pos, extfor 49 USE ions_base, ONLY : ions_vrescal, fricp, greasp, & 50 iforce, ndfrz, ions_shiftvar, ityp, & 51 atm, cdm, cdms, ions_cofmsub 52 USE cell_base, ONLY : at, bg, ainv, frich, & 53 greash, tpiba2, omega, alat, ibrav, & 54 celldm, h, hold, hnew, velh, & 55 wmass, press, iforceh, cell_force 56 USE local_pseudo, ONLY : allocate_local_pseudo 57 USE io_global, ONLY : stdout, ionode, ionode_id 58 USE dener, ONLY : detot 59 USE constants, ONLY : pi, k_boltzmann_au, au_ps 60 USE io_files, ONLY : psfile, pseudo_dir 61 USE wave_base, ONLY : wave_steepest, wave_verlet 62 USE wave_base, ONLY : wave_speed2, frice, grease 63 USE control_flags, ONLY : conv_elec, tconvthrs 64 USE check_stop, ONLY : check_stop_now 65 USE efcalc, ONLY : clear_nbeg, ef_force 66 USE ions_base, ONLY : zv, ions_vel 67 USE cp_electronic_mass, ONLY : emass, emass_cutoff, emass_precond 68 USE ions_positions, ONLY : tau0, taum, taup, taus, tausm, tausp, & 69 vels, velsm, velsp, ions_hmove, & 70 ions_move, fion, fionm 71 USE ions_nose, ONLY : gkbt, kbt, qnp, ndega, nhpcl, nhpdim, & 72 nhpbeg, nhpend, & 73 vnhp, xnhp0, xnhpm, xnhpp, & 74 atm2nhp, ions_nosevel, ions_noseupd, & 75 tempw, ions_nose_nrg, gkbt2nhp, & 76 ekin2nhp, anum2nhp 77 USE electrons_nose, ONLY : qne, ekincw, xnhe0, xnhep, xnhem, & 78 vnhe, electrons_nose_nrg, & 79 electrons_nose_shiftvar, & 80 electrons_nosevel, electrons_noseupd 81 USE pres_ai_mod, ONLY : P_ext, P_in, P_fin, pvar, volclu, & 82 surfclu, Surf_t, abivol, abisur 83 USE wavefunctions, ONLY : c0_bgrp, cm_bgrp, phi_bgrp 84 USE wannier_module, ONLY : allocate_wannier 85 USE cp_interfaces, ONLY : printout_new, move_electrons, newinit 86 USE cell_nose, ONLY : xnhh0, xnhhm, xnhhp, vnhh, temph, & 87 qnh, cell_nosevel, cell_noseupd, & 88 cell_nose_nrg, cell_nose_shiftvar 89 USE cell_base, ONLY : cell_kinene, cell_gamma, & 90 cell_move, cell_hmove 91 USE gvecw, ONLY : ecutwfc 92 USE gvect, ONLY : ecutrho 93 USE time_step, ONLY : delt, tps, dt2, twodelt 94 USE cp_interfaces, ONLY : cp_print_rho, nlfh, prefor, dotcsc 95 USE cp_main_variables, ONLY : acc, lambda, lambdam, lambdap, & 96 ema0bg, sfac, eigr, iprint_stdout, & 97 irb, taub, eigrb, rhog, rhos, & 98 rhor, bephi, becp_bgrp, nfi, idesc, & 99 drhor, drhog, bec_bgrp, dbec 100 USE autopilot, ONLY : event_step, event_index, & 101 max_event_step, restart_p 102 USE cell_base, ONLY : s_to_r, r_to_s 103 USE wannier_subroutines, ONLY : wannier_startup, wf_closing_options, & 104 ef_enthalpy 105 USE cp_interfaces, ONLY : writefile, eigs, strucf, phfacs 106 USE cp_interfaces, ONLY : ortho, elec_fakekine, calbec_bgrp, calbec, caldbec_bgrp 107 USE constraints_module, ONLY : check_constraint, remove_constr_force 108 USE cp_autopilot, ONLY : pilot 109 USE ions_nose, ONLY : ions_nose_allocate, ions_nose_shiftvar 110 USE orthogonalize_base, ONLY : updatc 111 USE control_flags, ONLY : force_pairing, tprint 112 USE mp, ONLY : mp_bcast, mp_sum 113 USE mp_global, ONLY : root_bgrp, intra_bgrp_comm, & 114 me_bgrp, inter_bgrp_comm, nbgrp, me_image 115 USE ldaU_cp, ONLY : lda_plus_u, vupsi 116 USE fft_base, ONLY : dfftp, dffts 117 USE london_module, ONLY : energy_london, force_london, stres_london 118 USE input_parameters, ONLY : tcpbo 119 USE funct, ONLY : dft_is_hybrid, start_exx, exx_is_active 120 USE funct, ONLY : dft_is_meta 121 ! 122 IMPLICIT NONE 123 ! 124 include 'laxlib.fh' 125 ! 126 ! ... input/output variables 127 ! 128 REAL(DP), INTENT(OUT) :: tau_out(3,nat) 129 REAL(DP), INTENT(OUT) :: fion_out(3,nat) 130 REAL(DP), INTENT(OUT) :: etot_out 131 ! 132 ! ... control variables 133 ! 134 LOGICAL :: tstop, tconv 135 LOGICAL :: tfile, tstdout 136 ! logical variable used to control printout 137 ! 138 ! ... forces on ions 139 ! 140 REAL(DP) :: maxfion, fion_tot(3) 141 ! 142 ! ... work variables 143 ! 144 REAL(DP) :: tempp, savee, saveh, savep, epot, epre, & 145 enow, econs, econt, fccc, ccc, bigr, dt2bye 146 REAL(DP) :: ekinc0, ekinp, ekinpr, ekinc 147 REAL(DP) :: temps(nsp) 148 REAL(DP) :: ekinh, temphc 149 REAL(DP) :: delta_etot 150 REAL(DP) :: ftmp, enb, enbi 151 INTEGER :: is, nacc, ia, j, iter, i, isa, ipos, iat, CYCLE_NOSE 152 INTEGER :: k, ii, l, m, iss 153 REAL(DP) :: hgamma(3,3), temphh(3,3) 154 REAL(DP) :: fcell(3,3) 155 REAL(DP) :: deltaP, ekincf 156 REAL(DP) :: stress_gpa(3,3), thstress(3,3), stress(3,3) 157 ! 158 REAL(DP), ALLOCATABLE :: usrt_tau0(:,:), usrt_taup(:,:), usrt_fion(:,:) 159 ! temporary array used to store unsorted positions and forces for 160 ! constrained dynamics 161 CHARACTER(LEN=3) :: labelw( nat ) 162 ! for force_pairing 163 INTEGER :: nspin_sub , i1, i2 164 ! pmass contains masses in atomic Hartree units 165 REAL(DP), ALLOCATABLE :: pmass(:) 166 REAL(DP), ALLOCATABLE :: forceh(:,:) 167 ! 168 REAL(DP) :: exx_start_thr 169 CALL start_clock( 'cpr_total' ) 170 ! 171 etot_out = 0.D0 172 enow = 1.D9 173 stress = 0.0D0 174 thstress = 0.0D0 175 ! moved to control_flags.f90 (Modules) 176 ! tfirst = .TRUE. 177 ! tlast = .FALSE. 178 nacc = 5 179 ! 180 if (dft_is_meta()) then 181 !HK/MCA : for SCAN0 calculation the initial SCAN has to converge better than the PBE -> PBE0 case 182 exx_start_thr = 1.E+1_DP 183 else 184 exx_start_thr = 1.E+2_DP 185 end if ! dft_is_meta 186 ! 187 ALLOCATE ( pmass (nsp) ) 188 pmass(1:nsp) = amass(1:nsp) * amu_au 189 nspin_sub = nspin 190 IF( force_pairing ) nspin_sub = 1 191 ! 192 ! ... Check for restart_p from Autopilot Feature Suite 193 ! 194 IF ( restart_p ) THEN 195 ! 196 ! ... do not add past nfi 197 ! 198 nomore = nomore 199 ! 200 END IF 201 ! 202 IF ( lda_plus_u ) ALLOCATE( forceh( 3, nat ) ) 203 ! 204 ! 205 !====================================================================== 206 ! 207 ! basic loop for molecular dynamics starts here 208 ! 209 !====================================================================== 210 ! 211 main_loop: DO 212 ! 213 CALL start_clock( 'main_loop' ) 214 CALL start_clock( 'cpr_md' ) 215 ! 216 dt2bye = dt2 / emass 217 nfi = nfi + 1 218 tlast = ( nfi == nomore ) .OR. tlast 219 tprint = ( MOD( nfi, iprint ) == 0 ) .OR. tlast !this can be set to .true. also by cp_autopilot in 'call pilot(nfi)', to compute velocities of the wfc in the last step of CG 220 tfile = ( MOD( nfi, iprint ) == 0 ) 221 tstdout = ( MOD( nfi, iprint_stdout ) == 0 ) .OR. tlast 222 ! 223 224 IF ( abivol ) THEN 225 IF ( pvar ) THEN 226 IF ( nfi .EQ. 1 ) THEN 227 deltaP = (P_fin - P_in) / DBLE(nomore) 228 P_ext = P_in 229 ELSE 230 P_ext = P_ext + deltaP 231 END IF 232 END IF 233 END IF 234 ! 235 IF ( ionode .AND. tstdout ) & 236 WRITE( stdout, '(/," * Physical Quantities at step:",I6)' ) nfi 237 ! 238 IF ( tnosee ) THEN 239 fccc = 1.D0 / ( 1.D0 + 0.5D0 * delt * vnhe ) 240 ELSE IF ( tsde ) THEN 241 fccc = 1.D0 242 ELSE 243 fccc = 1.D0 / ( 1.D0 + frice ) 244 END IF 245 ! 246 ! ... calculation of velocity of nose-hoover variables 247 ! 248 IF ( tnosep ) THEN 249 ! 250 CALL ions_nosevel( vnhp, xnhp0, xnhpm, delt, nhpcl, nhpdim ) 251 ! 252 END IF 253 ! 254 IF ( tnosee ) THEN 255 ! 256 CALL electrons_nosevel( vnhe, xnhe0, xnhem, delt ) 257 ! 258 END IF 259 ! 260 IF ( tnoseh ) THEN 261 ! 262 CALL cell_nosevel( vnhh, xnhh0, xnhhm, delt ) 263 ! 264 velh(:,:) = 2.D0 * ( h(:,:) - hold(:,:) ) / delt - velh(:,:) 265 ! 266 END IF 267 ! 268 IF ( (okvan .or. nlcc_any ) .AND. (tfor .OR. thdyn .OR. tfirst) ) THEN 269 ! 270 CALL initbox( tau0, alat, at, ainv, taub, irb ) 271 ! 272 CALL phbox( taub, iverbosity, eigrb ) 273 ! 274 END IF 275 ! 276 IF ( tfor .OR. thdyn ) THEN 277 ! 278 CALL phfacs( eigts1,eigts2,eigts3, eigr, mill, taus, dfftp%nr1,dfftp%nr2,dfftp%nr3, nat ) 279 ! 280 ! ... strucf calculates the structure factor sfac 281 ! 282 CALL strucf( sfac, eigts1, eigts2, eigts3, mill, dffts%ngm ) 283 ! 284 END IF 285 ! 286 IF ( thdyn ) THEN 287 ! 288 CALL formf( tfirst, eself ) 289 ! 290 END IF 291 ! 292 ! ... why this call ??? from Paolo Umari 293 ! 294 IF ( tefield .or. tefield2 ) THEN 295 ! 296 CALL calbec( 1, nsp, eigr, c0_bgrp, bec_bgrp ) ! ATTENZIONE 297 ! 298 END IF 299 ! 300 ! Autopilot (Dynamic Rules) Implimentation 301 ! 302 call pilot(nfi) 303 ! 304 IF ( ( tfor .OR. tfirst ) .AND. tefield ) CALL efield_update( eigr ) 305 IF ( ( tfor .OR. tfirst ) .AND. tefield2 ) CALL efield_update2( eigr ) 306 ! 307 ! ... pass ions information to plugins 308 ! 309 CALL plugin_init_ions( tau0 ) 310 ! 311 IF ( lda_plus_u ) then 312 ! forceh ! Forces on ions due to Hubbard U 313 forceh=0.0d0 314 ! vupsi ! potentials on electrons due to Hubbard U 315 vupsi=(0.0d0,0.0d0) 316 CALL new_ns(c0_bgrp,eigr,vkb,vupsi,forceh) 317 if ( mod(nfi,iprint).eq.0 ) call write_ns 318 endif 319 ! 320 !======================================================================= 321 ! 322 ! electronic degrees of freedom are updated here 323 ! 324 !======================================================================= 325 ! 326 IF( force_pairing ) THEN 327 c0_bgrp(:,iupdwn(2):nbsp) = c0_bgrp(:,1:nupdwn(2)) 328 cm_bgrp(:,iupdwn(2):nbsp) = cm_bgrp(:,1:nupdwn(2)) 329 phi_bgrp(:,iupdwn(2):nbsp) = phi_bgrp(:,1:nupdwn(2)) 330 lambda(:,:, 2) = lambda(:,:, 1) 331 ENDIF 332 ! 333 ! ... fake electronic kinetic energy 334 ! 335 IF ( .NOT. tcg ) THEN 336 ! 337 ekincf = 0.0d0 338 339 CALL elec_fakekine( ekincf, ema0bg, emass, cm_bgrp, c0_bgrp, ngw, nbsp_bgrp, 1, delt ) 340 ! 341 END IF 342 ! 343 CALL move_electrons( nfi, tfirst, tlast, bg(:,1), bg(:,2), bg(:,3), & 344 fion, c0_bgrp, cm_bgrp, phi_bgrp, & 345 enthal, enb, enbi, fccc, ccc, dt2bye, stress, .false. ) 346 ! 347 IF (lda_plus_u) fion = fion + forceh 348 ! 349 ! DFT+D (Grimme) dispersion energy, forces (factor 0.5 converts to Ha/a.u.) 350 ! 351 IF ( llondon ) THEN 352 ALLOCATE( usrt_tau0( 3, nat )) 353 usrt_tau0(:,:) = tau0(:,:)/alat 354 delta_etot = 0.5_dp*energy_london (alat, nat,ityp,at,bg, usrt_tau0) 355 etot = etot + delta_etot 356 enthal=enthal+delta_etot 357 IF ( tfor ) THEN 358 ALLOCATE( usrt_fion( 3, nat ) ) 359 usrt_fion = 0.5_dp*force_london ( alat, nat,ityp, at,bg, usrt_tau0 ) 360 fion(:,:) = fion(:,:) + usrt_fion(:,:) 361 DEALLOCATE (usrt_fion) 362 END IF 363 IF ( tpre ) stress = stress + 0.5_dp * stres_london ( alat , nat , & 364 ityp , at , bg , usrt_tau0 , omega ) 365 DEALLOCATE ( usrt_tau0 ) 366 END IF 367 ! 368 IF ( tpre ) THEN 369 ! 370 CALL nlfh( stress, bec_bgrp, dbec, lambda, idesc ) 371 ! 372 CALL ions_thermal_stress( stress, thstress, pmass, omega, h, vels, nat, ityp ) 373 ! 374 IF (tstdout) THEN 375 WRITE(stdout,'(5X,"Pressure of Nuclei (GPa)",F20.5,I7)') & 376 (thstress(1,1)+thstress(2,2)+thstress(3,3))/3.0_DP * au_gpa, nfi 377 WRITE(stdout,'(5X,"Pressure Total (GPa)",F20.5,I7)') & 378 (stress(1,1)+stress(2,2)+stress(3,3))/3.0_DP * au_gpa , nfi 379 END IF 380 ! 381 END IF 382 ! 383 ! 384 !======================================================================= 385 ! 386 ! verlet algorithm 387 ! 388 ! loop which updates cell parameters and ionic degrees of freedom 389 ! hnew=h(t+dt) is obtained from hold=h(t-dt) and h=h(t) 390 ! tausp=pos(t+dt) from tausm=pos(t-dt) taus=pos(t) h=h(t) 391 ! 392 ! guessed displacement of ions 393 !======================================================================= 394 ! 395 hgamma(:,:) = 0.D0 396 ! 397 IF ( thdyn ) THEN 398 ! 399 CALL cell_force( fcell, ainv, stress, omega, press, wmass ) 400 ! 401 CALL cell_move( hnew, h, hold, delt, iforceh, & 402 fcell, frich, tnoseh, vnhh, velh, tsdc ) 403 ! 404 velh(:,:) = ( hnew(:,:) - hold(:,:) ) / twodelt 405 ! 406 CALL cell_gamma( hgamma, ainv, h, velh ) 407 ! 408 END IF 409 ! 410 ! BS : Initialization of additional cycles for the Nose thermostat ... 411 ! 412 IF (tnosep) CYCLE_NOSE=0 413 ! 414444 IF ( tfor ) THEN 415 ! 416 IF ( lwf ) CALL ef_force( fion, ityp, nat, zv ) 417 ! 418 IF( textfor ) THEN 419 ! 420 FORALL( ia = 1:nat ) fion(:,ia) = fion(:,ia) + extfor(:,ia) 421 ! 422 fion_tot(:) = SUM( fion(:,:), DIM = 2 ) / DBLE( nat ) 423 ! 424 FORALL( ia = 1:nat ) fion(:,ia) = fion(:,ia) - fion_tot(:) 425 ! 426 END IF 427 ! 428 IF ( remove_rigid_rot ) & 429 CALL remove_tot_torque( nat, tau0, pmass(ityp(:)), fion ) 430 ! 431 IF ( lconstrain ) THEN 432 ! 433 IF ( ionode ) THEN 434 ! 435 ALLOCATE( usrt_tau0( 3, nat ) ) 436 ALLOCATE( usrt_taup( 3, nat ) ) 437 ALLOCATE( usrt_fion( 3, nat ) ) 438 ! 439 usrt_tau0(:,:) = tau0(:,:) 440 usrt_fion(:,:) = fion(:,:) 441 ! 442 ! ... we first remove the component of the force along the 443 ! ... constrain gradient (this constitutes the initial guess 444 ! ... for the lagrange multiplier) 445 ! 446 CALL remove_constr_force( nat, usrt_tau0, if_pos, ityp, 1.D0, usrt_fion ) 447 ! 448 fion(:,:) = usrt_fion(:,:) 449 ! 450 END IF 451 ! 452 CALL mp_bcast( fion, ionode_id, intra_bgrp_comm ) 453 ! 454 END IF 455 ! 456 ! 457 ! ... call void routine for user define/ plugin patches on external forces 458 ! 459 CALL plugin_ext_forces() 460 ! 461 ! 462 CALL ions_move( tausp, taus, tausm, iforce, pmass, fion, ainv, & 463 delt, ityp, nat, fricp, hgamma, vels, tsdp, tnosep, & 464 fionm, vnhp, velsp, velsm, nhpcl, nhpdim, atm2nhp ) 465 ! 466 IF ( lconstrain ) THEN 467 ! 468 ! ... constraints are imposed here 469 ! 470 IF ( ionode ) THEN 471 ! 472 CALL s_to_r( tausp, taup, nat, hnew ) 473 ! 474 usrt_taup(:,:) = taup(:,:) 475 ! 476 CALL check_constraint( nat, usrt_taup, usrt_tau0, usrt_fion, & 477 if_pos, ityp, 1.D0, delt, amu_au ) 478 ! 479 taup(:,:) = usrt_taup(:,:) 480 fion(:,:) = usrt_fion(:,:) 481 ! 482 DEALLOCATE( usrt_tau0, usrt_taup, usrt_fion ) 483 ! 484 END IF 485 ! 486 CALL mp_bcast( taup, ionode_id, intra_bgrp_comm ) 487 CALL mp_bcast( fion, ionode_id, intra_bgrp_comm ) 488 ! 489 CALL r_to_s( taup, tausp, nat, ainv ) 490 ! 491 END IF 492 ! 493 CALL ions_cofmass( tausp, pmass, nat, ityp, cdm ) 494 ! 495 IF ( ndfrz == 0 ) & 496 CALL ions_cofmsub( tausp, iforce, nat, cdm, cdms ) 497 ! 498 CALL s_to_r( tausp, taup, nat, hnew ) 499 ! 500 END IF 501 ! 502 !-------------------------------------------------------------------------- 503 ! initialization with guessed positions of ions 504 !-------------------------------------------------------------------------- 505 ! 506 ! ... if thdyn=true g vectors and pseudopotentials are recalculated for 507 ! ... the new cell parameters 508 ! 509 ! 510 IF ( tfor .OR. thdyn ) THEN 511 ! 512 IF ( .NOT.tnosep .OR. CYCLE_NOSE.EQ.0 ) THEN 513 ! 514 IF ( thdyn ) THEN 515 ! 516 hold = h 517 h = hnew 518 ! 519 IF( nbgrp > 1 ) THEN 520 CALL mp_bcast( h, 0, inter_bgrp_comm ) 521 END IF 522 ! 523 CALL newinit( h, iverbosity ) 524 ! 525 CALL newnlinit() 526 ! 527 ELSE 528 ! 529 hold = h 530 ! 531 END IF 532 ! 533 END IF 534 ! 535 ! ... phfac calculates eigr 536 ! 537 CALL phfacs( eigts1,eigts2,eigts3, eigr, mill, tausp, dfftp%nr1,dfftp%nr2,dfftp%nr3, nat ) 538 ! ... prefor calculates vkb 539 ! 540 CALL prefor( eigr, vkb ) 541 ! 542 END IF 543 ! 544 !-------------------------------------------------------------------------- 545 ! imposing the orthogonality 546 !-------------------------------------------------------------------------- 547 ! 548 ! In case of tnosep = .true., the orthonormality is done only with the most updated 549 ! atomic coordinates coming out of the CYCLE_NOSE loop 550 ! 551 IF ( .NOT.tnosep .OR. CYCLE_NOSE.EQ.2 ) THEN 552 ! 553 IF ( .NOT. tcg ) THEN 554 ! 555 IF ( tortho ) THEN 556 ! 557 CALL ortho( eigr, cm_bgrp, phi_bgrp, lambda, idesc, bigr, iter, ccc, bephi, becp_bgrp ) 558 ! 559 ELSE 560 ! 561 CALL gram_bgrp( vkb, bec_bgrp, nkb, cm_bgrp, ngw ) 562 ! 563 IF ( iverbosity > 2 ) CALL dotcsc( eigr, cm_bgrp, ngw, nbsp_bgrp ) 564 ! 565 END IF 566 ! 567 ! correction to displacement of ions 568 ! 569 IF ( iverbosity > 1 ) CALL laxlib_print_matrix( lambda, idesc, nbsp, 9, nudx, 1.D0, ionode, stdout ) 570 ! 571 IF ( tortho ) THEN 572 CALL updatc( ccc, lambda, phi_bgrp, bephi, becp_bgrp, bec_bgrp, cm_bgrp, idesc ) 573 END IF 574 ! 575 IF( force_pairing ) THEN 576 c0_bgrp(:,iupdwn(2):nbsp) = c0_bgrp(:,1:nupdwn(2)) 577 cm_bgrp(:,iupdwn(2):nbsp) = cm_bgrp(:,1:nupdwn(2)) 578 phi_bgrp(:,iupdwn(2):nbsp) = phi_bgrp(:,1:nupdwn(2)) 579 lambda(:,:, 2) = lambda(:,:, 1) 580 ENDIF 581 ! 582 CALL calbec_bgrp( 1, nsp, eigr, cm_bgrp, bec_bgrp, 1 ) 583 ! 584 IF ( tpre ) THEN 585 CALL caldbec_bgrp( eigr, cm_bgrp, dbec, idesc ) 586 END IF 587 ! 588 IF ( iverbosity > 1 ) CALL dotcsc( eigr, cm_bgrp, ngw, nbsp_bgrp ) 589 ! 590 END IF 591 ! 592 END IF !(.NOT.tnosep.OR.CYCLE_NOSE.EQ.2) 593 ! 594 !-------------------------------------------------------------------------- 595 ! temperature monitored and controlled 596 !-------------------------------------------------------------------------- 597 ! 598 ekinp = 0.D0 599 ekinpr = 0.D0 600 tempp = 0.D0 601 temps = 0.D0 602 ekinc0 = 0.0d0 603 ekinc = 0.0d0 604 ! 605 ! 606 ! ... ionic kinetic energy and temperature 607 ! 608 IF ( tfor ) THEN 609 ! 610 CALL ions_vel( vels, tausp, tausm, delt ) 611 ! 612 CALL ions_kinene( ekinp, vels, nat, ityp, hold, pmass ) 613 ! 614 CALL ions_temp( tempp, temps, ekinpr, vels, nsp, na, nat, ityp, & 615 hold, pmass, ndega, nhpdim, atm2nhp, ekin2nhp ) 616 ! 617 END IF 618 ! 619 ! ... fake electronic kinetic energy 620 ! 621 IF ( .NOT. tcg ) THEN 622 ! 623 CALL elec_fakekine( ekinc0, ema0bg, emass, c0_bgrp, cm_bgrp, ngw, nbsp_bgrp, 1, delt ) 624 ! 625 ekinc0 = (ekinc0 + ekincf)*0.5d0 626 ! 627 ekinc = ekinc0 628 ! 629 END IF 630 ! 631 ! ... fake cell-parameters kinetic energy 632 ! 633 ekinh = 0.D0 634 ! 635 IF ( thdyn ) THEN 636 ! 637 CALL cell_kinene( ekinh, temphh, velh ) 638 ! 639 END IF 640 ! 641 IF ( COUNT( iforceh == 1 ) > 0 ) THEN 642 ! 643 temphc = 2.D0 / k_boltzmann_au * ekinh / DBLE( COUNT( iforceh == 1 ) ) 644 ! 645 ELSE 646 ! 647 temphc = 0.D0 648 ! 649 END IF 650 ! 651 ! ... udating nose-hoover friction variables 652 ! 653 IF ( tnosep ) CALL ions_noseupd( xnhpp, xnhp0, xnhpm, delt, qnp, & 654 ekin2nhp, gkbt2nhp, vnhp, kbt, & 655 nhpcl, nhpdim, nhpbeg, nhpend ) 656 ! 657 IF ( tnosee ) CALL electrons_noseupd( xnhep, xnhe0, xnhem, & 658 delt, qne, ekinc, ekincw, vnhe ) 659 ! 660 IF ( tnoseh ) CALL cell_noseupd( xnhhp, xnhh0, xnhhm, & 661 delt, qnh, temphh, temph, vnhh ) 662 ! 663 !================================================================= 664 ! BS : Additional cycles for the Nose thermostat ... 665 IF(tnosep) CYCLE_NOSE=CYCLE_NOSE+1 666 IF(tnosep .AND. (CYCLE_NOSE .LE. 2) ) GO TO 444 667 !================================================================= 668 ! 669 ! ... warning: thdyn and tcp/tcap are not compatible yet!!! 670 ! 671 IF ( tcp .AND. tfor .AND. .NOT.thdyn ) THEN 672 ! 673 IF ( tempp > (tempw+tolp) .OR. & 674 tempp < (tempw-tolp) .AND. tempp /= 0.D0 ) THEN 675 ! 676 CALL ions_vrescal( .false., tempw, tempp, taup, & 677 tau0, taum, nat, ityp, fion, iforce, pmass, delt ) 678 CALL r_to_s( taup, tausp, nat, ainv ) 679 ! 680 END IF 681 ! 682 END IF 683 ! 684 IF ( tprint ) THEN 685 ! 686 IF( tortho ) THEN 687 ! 688 IF( force_pairing ) THEN 689 lambda(:, :, 2) = lambda(:, :, 1) 690 lambdap(:, :, 2) = lambdap(:, :, 1) 691 WRITE( stdout, '("Occupations in CPR:")' ) 692 WRITE( stdout, '(10F9.6)' ) ( f(i), i = 1, nbspx ) 693 END IF 694 ! 695 CALL eigs( nfi, lambdap, lambda, idesc ) 696 ! 697 ELSE 698 ! 699 WRITE( stdout, '("NOTE: eigenvalues are not computed without ortho")' ) 700 ! 701 END IF 702 ! 703 END IF 704 ! 705 IF ( lwf ) CALL ef_enthalpy( enthal, tau0 ) 706 ! 707 IF ( tens .AND. tprint ) THEN 708 ! 709 WRITE( stdout, '("Occupations :")' ) 710 WRITE( stdout, '(10F9.6)' ) ( f(i), i = 1, nbsp ) 711 ! 712 END IF 713 ! 714 epot = eht + epseu + exc 715 ! 716 IF ( .NOT. tcg ) THEN 717 ! 718 econs = ekinp + ekinh + enthal 719 econt = econs + ekinc 720 ! 721 ELSE 722 ! 723 IF ( .NOT. tens ) THEN 724 ! 725 econs = ekinp + etot 726 atot = etot 727 econt = econs 728 ! 729 ELSE 730 ! 731 gibbsfe = atot 732 econs = ekinp + atot 733 econt = econs 734 ! 735 END IF 736 ! 737 END IF 738 ! 739 ! ... add energies of thermostats 740 ! 741 IF ( tnosep ) & 742 econt = econt + ions_nose_nrg( xnhp0, vnhp, qnp, & 743 gkbt2nhp, kbt, nhpcl, nhpdim ) 744 IF ( tnosee ) & 745 econt = econt + electrons_nose_nrg( xnhe0, vnhe, qne, ekincw ) 746 IF ( tnoseh ) & 747 econt = econt + cell_nose_nrg( qnh, xnhh0, vnhh, temph, iforceh ) 748 ! 749 tps = tps + delt * au_ps 750 ! 751 if (abivol) etot = etot - P_ext*volclu 752 if (abisur) etot = etot - Surf_t*surfclu 753 ! 754 IF ( tstdout) CALL spinsq ( c0_bgrp, bec_bgrp, rhor ) 755 ! 756 CALL printout_new( nfi, tfirst, tfile, tprint, tps, hold, stress, & 757 tau0, vels, fion, ekinc, temphc, tempp, temps, etot, & 758 enthal, econs, econt, vnhh, xnhh0, vnhp, xnhp0, vnhe, xnhe0, atot, & 759 ekin, epot, tprnfor, tpre, tstdout ) 760 ! 761 if (abivol) etot = etot + P_ext*volclu 762 if (abisur) etot = etot + Surf_t*surfclu 763 ! 764 IF( tfor ) THEN 765 ! 766 ! ... new variables for next step 767 ! 768 CALL ions_shiftvar( taup, tau0, taum ) ! real positions 769 CALL ions_shiftvar( tausp, taus, tausm ) ! scaled positions 770 CALL ions_shiftvar( velsp, vels, velsm ) ! scaled velocities 771 ! 772 IF ( tnosep ) CALL ions_nose_shiftvar( xnhpp, xnhp0, xnhpm ) 773 IF ( tnosee ) CALL electrons_nose_shiftvar( xnhep, xnhe0, xnhem ) 774 IF ( tnoseh ) CALL cell_nose_shiftvar( xnhhp, xnhh0, xnhhm ) 775 ! 776 IF( nbgrp > 1 ) THEN 777 CALL mp_bcast( tau0, 0, inter_bgrp_comm ) 778 CALL mp_bcast( taus, 0, inter_bgrp_comm ) 779 CALL mp_bcast( vels, 0, inter_bgrp_comm ) 780 CALL mp_bcast( xnhp0, 0, inter_bgrp_comm ) 781 CALL mp_bcast( xnhe0, 0, inter_bgrp_comm ) 782 CALL mp_bcast( xnhh0, 0, inter_bgrp_comm ) 783 END IF 784 ! 785 END IF 786 ! 787 ekincm = ekinc0 788 ! 789 ! ... cm=c(t+dt) c0=c(t) 790 ! 791 IF( .NOT. tcg ) THEN 792 ! 793 CALL dswap( 2*SIZE( c0_bgrp ), c0_bgrp, 1, cm_bgrp, 1 ) 794 ! 795 ELSE 796 ! 797 CALL cg_update( tfirst, nfi, c0_bgrp ) 798 ! 799 IF ( tfor .AND. .NOT. tens .AND. tprint ) THEN 800 ! 801 ! ... in this case optimize c0 and lambda for smooth 802 ! ... restart with CP 803 ! 804 IF ( okvan .or. nlcc_any ) THEN 805 CALL initbox( tau0, alat, at, ainv, taub, irb ) 806 CALL phbox( taub, iverbosity, eigrb ) 807 END IF 808 CALL r_to_s( tau0, taus, nat, ainv ) 809 CALL phfacs( eigts1,eigts2,eigts3, eigr, mill, taus, dfftp%nr1,dfftp%nr2,dfftp%nr3, nat ) 810 CALL strucf( sfac, eigts1, eigts2, eigts3, mill, dffts%ngm ) 811 ! 812 IF ( thdyn ) CALL formf( tfirst, eself ) 813 IF ( tefield ) CALL efield_update( eigr ) 814 IF ( tefield2 ) CALL efield_update2( eigr ) 815 ! 816 CALL plugin_init_ions( tau0 ) 817 ! 818 lambdam = lambda 819 ! 820 CALL move_electrons( nfi, tfirst, tlast, bg(:,1), bg(:,2), bg(:,3),& 821 fion, c0_bgrp, cm_bgrp, phi_bgrp, enthal, enb,& 822 enbi, fccc, ccc, dt2bye, stress,.true. ) 823 ! 824 END IF 825 ! 826 END IF 827 ! 828 ! ... now: cm=c(t) c0=c(t+dt) 829 ! ... and, if tcg == .true. : 830 ! ... c0old=c(t),c0=c(t+dt) 831 ! 832 tfirst = .FALSE. 833 ! 834 CALL stop_clock( 'main_loop' ) 835 ! 836 ! ... write on file ndw each isave 837 ! 838 IF ( ( MOD( nfi, isave ) == 0 ) .AND. ( nfi < nomore ) ) THEN 839 ! 840 IF ( tcg ) THEN 841 ! 842 CALL writefile( h, hold ,nfi, c0_bgrp, c0old, taus, tausm, & 843 vels, velsm, acc, lambda, lambdam, idesc, xnhe0, xnhem, & 844 vnhe, xnhp0, xnhpm, vnhp, nhpcl,nhpdim,ekincm, xnhh0,& 845 xnhhm, vnhh, velh, fion, tps, z0t, f, rhor ) 846 ! 847 ELSE 848 ! 849 CALL writefile( h, hold, nfi, c0_bgrp, cm_bgrp, taus, & 850 tausm, vels, velsm, acc, lambda, lambdam, idesc, xnhe0, & 851 xnhem, vnhe, xnhp0, xnhpm, vnhp, nhpcl, nhpdim, ekincm,& 852 xnhh0, xnhhm, vnhh, velh, fion, tps, z0t, f, rhor ) 853 ! 854 END IF 855 ! 856 END IF 857 ! 858 epre = enow 859 enow = etot 860 ! 861 frice = frice * grease 862 fricp = fricp * greasp 863 frich = frich * greash 864 ! 865 !====================================================================== 866 ! 867 delta_etot = ABS( epre - enow ) 868 ! 869 !exx_wf related 870 !The following criteria is used to turn on exact exchange calculation when 871 !GGA energy is converged up to 100 times of the input etot convergence thereshold 872 ! 873 IF( .NOT.exx_is_active().AND.dft_is_hybrid().AND.tconvthrs%active ) THEN 874 ! 875 IF(delta_etot.LT.tconvthrs%derho*exx_start_thr) THEN 876 ! 877 WRITE(stdout,'(/,3X,"Exact Exchange is turned on ...")') 878 ! 879 CALL start_exx() 880 ! 881 END IF 882 ! 883 END IF 884 ! 885 tstop = check_stop_now() .OR. tlast 886 ! 887 tconv = .FALSE. 888 ! 889 IF ( tconvthrs%active ) THEN 890 ! 891 ! ... electrons 892 ! 893 tconv = ( ekinc < tconvthrs%ekin .AND. delta_etot < tconvthrs%derho ) 894 ! 895 IF ( tfor ) THEN 896 ! 897 ! ... ions 898 ! 899 maxfion = MAXVAL( ABS( fion(:,1:nat) ) ) 900 ! 901 tconv = tconv .AND. ( maxfion < tconvthrs%force ) 902 ! 903 END IF 904 ! 905 END IF 906 ! 907 ! ... in the case cp-wf the check on convergence is done starting 908 ! ... from the second step 909 ! 910 IF ( lwf .AND. tfirst ) tconv = .FALSE. 911 ! 912 IF ( tconv ) THEN 913 ! 914 tlast = .TRUE. 915 ! 916 IF ( ionode ) THEN 917 ! 918 WRITE( stdout, & 919 & "(/,3X,'MAIN:',10X,'EKINC (thr)', & 920 & 10X,'DETOT (thr)',7X,'MAXFORCE (thr)')" ) 921 WRITE( stdout, "(3X,'MAIN: ',3(D14.6,1X,D8.1))" ) & 922 ekinc, tconvthrs%ekin, delta_etot, & 923 tconvthrs%derho, 0.D0, tconvthrs%force 924 WRITE( stdout, & 925 "(3X,'MAIN: convergence achieved for system relaxation')" ) 926 ! 927 END IF 928 ! 929 END IF 930 ! 931 IF ( lwf ) & 932 CALL wf_closing_options( nfi, c0_bgrp, cm_bgrp, bec_bgrp, eigr, eigrb,& 933 taub, irb, ibrav, bg(:,1), bg(:,2), bg(:,3), & 934 taus, tausm, vels, & 935 velsm, acc, lambda, lambdam, idesc, xnhe0, xnhem, & 936 vnhe, xnhp0, xnhpm, vnhp, nhpcl, nhpdim, & 937 ekincm, xnhh0, xnhhm, vnhh, velh, ecutrho, & 938 ecutwfc,delt,celldm, fion, tps, z0t, f, rhor ) 939 ! 940 IF ( tstop ) EXIT main_loop 941 ! 942 CALL stop_clock( 'cpr_md' ) 943 ! 944 END DO main_loop 945 ! 946 !===================== end of main loop of molecular dynamics =============== 947 ! 948 DEALLOCATE ( pmass ) 949 ! ... Here copy relevant physical quantities into the output arrays/variables 950 ! 951 etot_out = etot 952 ! 953 tau_out(:,:) = tau0(:,:) 954 fion_out(:,:) = fion(:,:) 955 ! 956 conv_elec = .TRUE. 957 ! 958 IF ( tcg ) cm_bgrp = c0old 959 ! 960 CALL writefile( h, hold, nfi, c0_bgrp, cm_bgrp, taus, tausm, & 961 vels, velsm, acc, lambda, lambdam, idesc, xnhe0, xnhem, vnhe, & 962 xnhp0, xnhpm, vnhp, nhpcl,nhpdim,ekincm, xnhh0, xnhhm, & 963 vnhh, velh, fion, tps, z0t, f, rhor ) 964 ! 965 IF( iverbosity > 1 ) CALL laxlib_print_matrix( lambda, idesc, nbsp, nbsp, nudx, 1.D0, ionode, stdout ) 966 ! 967 IF (lda_plus_u) DEALLOCATE( forceh ) 968 969 ! 970 CALL stop_clock( 'cpr_total' ) ! BS 971 ! 972 RETURN 973 ! 974END SUBROUTINE cprmain 975! 976!---------------------------------------------------------------------------- 977SUBROUTINE terminate_run() 978 !---------------------------------------------------------------------------- 979 ! 980 USE io_global, ONLY : stdout, ionode 981 USE control_flags, ONLY : ts_vdw, thdyn, tortho 982 USE cg_module, ONLY : tcg, print_clock_tcg 983 USE ldaU_cp, ONLY : lda_plus_u 984 USE mp, ONLY : mp_report 985 USE control_flags, ONLY : lwf, lwfpbe0nscf 986 USE tsvdw_module, ONLY : tsvdw_finalize 987 USE exx_module, ONLY : exx_finalize 988 USE funct, ONLY : dft_is_hybrid, exx_is_active 989 ! 990 IMPLICIT NONE 991 ! 992 ! ... print statistics 993 ! 994 CALL printacc() 995 ! 996 !============================================================== 997 WRITE( stdout, '(/5x,"Called by MAIN_LOOP:")' ) 998 CALL print_clock( 'total_time' ) 999 CALL print_clock( 'initialize' ) 1000 CALL print_clock( 'main_loop' ) 1001 CALL print_clock( 'cpr_total' ) 1002 !============================================================== 1003 WRITE( stdout, '(/5x,"Called by INIT_RUN:")' ) 1004 CALL print_clock( 'init_readfile' ) 1005 IF( lwf ) CALL print_clock( 'wf_start' ) 1006 !============================================================== 1007 WRITE( stdout, '(/5x,"Called by CPR:")' ) 1008 CALL print_clock( 'cpr_md' ) 1009 CALL print_clock( 'move_electrons' ) 1010 CALL print_clock( 'move_ion' ) 1011 IF( lwf ) CALL print_clock( 'wf_close_opt' ) ! wf_close_opt = wf_1 + wf_2 1012 !============================================================== 1013 IF( lwf ) THEN 1014 WRITE( stdout, '(/5x,"Called by WANNIER_MODULES:")' ) 1015 CALL print_clock( 'wf_start' ) 1016 CALL print_clock( 'wf_init' ) 1017 CALL print_clock( 'wf_close_opt' ) ! wf_close_opt = wf_1 + wf_2 1018 CALL print_clock( 'wf_1' ) 1019 CALL print_clock( 'wf_2' ) 1020 CALL print_clock( 'ddyn_u' ) 1021 CALL print_clock( 'ortho_u' ) 1022 END IF 1023 !============================================================== 1024 !exx_wf related 1025 IF ( dft_is_hybrid().AND.exx_is_active() ) THEN 1026 ! 1027 WRITE( stdout, '(/5x,"Called by EXACT_EXCHANGE:")' ) 1028 CALL print_clock('exact_exchange') ! total time for exx 1029 CALL print_clock('self_v') ! calculation self potential 1030 CALL print_clock('getpairv') ! calculation pair potential 1031 CALL print_clock('exx_gs_setup') ! calculation 1032 CALL print_clock('exx_pairs') ! calculation 1033 CALL print_clock('r_orbital') ! communication 1034 CALL print_clock('totalenergy') ! communication 1035 CALL print_clock('vl2vg') ! communication 1036 CALL print_clock('send_psi') ! communication 1037 CALL print_clock('sendv') ! communication 1038 CALL print_clock('send_psi_barrier') ! communication 1039 CALL print_clock('send_psi_wait') ! communication 1040 !CALL print_clock('sendv_barrier') ! communication 1041 CALL print_clock('getvofr') 1042 CALL print_clock('getvofr_qlm') 1043 CALL print_clock('getvofr_bound') 1044 CALL print_clock('getvofr_geterho') 1045 CALL print_clock('getvofr_hpotcg') 1046 !CALL print_clock('hpotcg') 1047 !CALL print_clock('getqlm') 1048 !CALL print_clock('exx_bound') 1049 !CALL print_clock('lapmv') 1050 CALL print_clock('exx_cell_derv') 1051 ! 1052 CALL exx_finalize() ! deallocate all arrays 1053 ! 1054 END IF 1055 !============================================================== 1056 IF (thdyn) THEN 1057 WRITE( stdout, '(/5x,"Called by CELL_DYNAMICS:")' ) 1058 CALL print_clock( 'formf' ) 1059 END IF 1060 1061 WRITE( stdout, '(/5x,"Called by move_electrons:")' ) 1062 CALL print_clock( 'rhoofr' ) 1063 CALL print_clock( 'vofrho' ) 1064 CALL print_clock( 'dforce' ) 1065 CALL print_clock( 'calphi' ) 1066 CALL print_clock( 'newd' ) 1067 CALL print_clock( 'nlfl' ) 1068 1069 IF (lda_plus_u) WRITE( stdout, '(/5x,"Called by new_ns:")' ) 1070 CALL print_clock( 'new_ns:forc' ) 1071 CALL print_clock( 'projwfc_hub' ) 1072 CALL print_clock( 'dndtau' ) 1073 1074 IF (tortho) WRITE( stdout, '(/5x,"Called by ortho:")' ) 1075 IF (tortho) THEN 1076 CALL print_clock( 'ortho_iter' ) 1077 CALL print_clock( 'rsg' ) 1078 CALL print_clock( 'rhoset' ) 1079 CALL print_clock( 'sigset' ) 1080 CALL print_clock( 'tauset' ) 1081 CALL print_clock( 'ortho' ) 1082 CALL print_clock( 'updatc' ) 1083 ELSE 1084 CALL print_clock( 'gram' ) 1085 END IF 1086 1087 WRITE( stdout, '(/5x,"Small boxes:")' ) 1088 CALL print_clock( 'rhov' ) 1089 CALL print_clock( 'fftb' ) 1090 CALL print_clock( 'set_cc' ) 1091 CALL print_clock( 'forcecc' ) 1092 1093 WRITE( stdout, '(/5x,"Low-level routines:")' ) 1094 CALL print_clock( 'prefor' ) 1095 CALL print_clock( 'nlfq' ) 1096 CALL print_clock( 'nlsm1' ) 1097 CALL print_clock( 'nlsm2' ) 1098 CALL print_clock( 'fft' ) 1099 CALL print_clock( 'ffts' ) 1100 CALL print_clock( 'fftw' ) 1101 CALL print_clock( 'fft_scatt_xy' ) 1102 CALL print_clock( 'fft_scatt_yz' ) 1103 CALL print_clock( 'fft_scatt_tg' ) 1104 CALL print_clock( 'betagx' ) 1105 CALL print_clock( 'qradx' ) 1106 CALL print_clock( 'tmp_clk1' ) 1107 CALL print_clock( 'tmp_clk2' ) 1108 CALL print_clock( 'tmp_clk3' ) 1109 CALL print_clock( 'gram' ) 1110 CALL print_clock( 'nlinit' ) 1111 CALL print_clock( 'init_dim' ) 1112 CALL print_clock( 'newnlinit' ) 1113 CALL print_clock( 'from_scratch' ) 1114 CALL print_clock( 'from_restart' ) 1115 CALL print_clock( 'new_ns' ) 1116 CALL print_clock( 'strucf' ) 1117 CALL print_clock( 'calbec' ) 1118!============================================================== 1119 IF (ts_vdw) THEN 1120 WRITE( stdout, '(/5x,"Called by tsvdw:")' ) 1121 CALL print_clock( 'ts_vdw' ) 1122 CALL print_clock( 'tsvdw_pair' ) 1123 CALL print_clock( 'tsvdw_rhotot' ) 1124 CALL print_clock( 'tsvdw_screen' ) 1125 CALL print_clock( 'tsvdw_veff' ) 1126 CALL print_clock( 'tsvdw_dveff' ) 1127 CALL print_clock( 'tsvdw_energy' ) 1128 CALL print_clock( 'tsvdw_wfforce' ) 1129 ! 1130 CALL tsvdw_finalize() 1131 END IF 1132 ! 1133 IF (tcg) call print_clock_tcg() 1134 ! 1135 CALL plugin_clock() 1136 ! 1137 CALL mp_report() 1138 ! 1139END SUBROUTINE terminate_run 1140