1! cp_autopilot.f90 2!******************************************************************************** 3! cp_autopilot.f90 Copyright (c) 2005 Targacept, Inc. 4!******************************************************************************** 5! The Autopilot Feature suite is a user level enhancement that enables the 6! following features: 7! automatic restart of a job; 8! preconfiguration of job parameters; 9! on-the-fly changes to job parameters; 10! and pausing of a running job. 11! 12! For more information, see AUTOPILOT in document directory. 13! 14! This program is free software; you can redistribute it and/or modify it under 15! the terms of the GNU General Public License as published by the Free Software 16! Foundation; either version 2 of the License, or (at your option) any later version. 17! This program is distributed in the hope that it will be useful, but WITHOUT ANY 18! WARRANTY; without even the implied warranty of MERCHANTABILITY FOR A PARTICULAR 19! PURPOSE. See the GNU General Public License at www.gnu.or/copyleft/gpl.txt for 20! more details. 21! 22! THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. 23! EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES 24! PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, 25! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 26! FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND THE 27! PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, 28! YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. 29! 30! IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING, 31! WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE 32! THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY 33! GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR 34! INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA 35! BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A 36! FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER 37! OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. 38! 39! You should have received a copy of the GNU General Public License along with 40! this program; if not, write to the 41! Free Software Foundation, Inc., 42! 51 Franklin Street, 43! Fifth Floor, 44! Boston, MA 02110-1301, USA. 45! 46! Targacept's address is 47! 200 East First Street, Suite 300 48! Winston-Salem, North Carolina USA 27101-4165 49! Attn: Molecular Design. 50! Email: atp@targacept.com 51! 52! This work was supported by the Advanced Technology Program of the 53! National Institute of Standards and Technology (NIST), Award No. 70NANB3H3065 54! 55!******************************************************************************** 56 57!******************************************************************************** 58! 04/2018 59! Implemented switch from CG to Verlet and from Verlet to CG. 60! The code automagically set the correct wfc velocity after the last CG step, 61! and inizialize again the mass of the electron mu(k) after the CG. 62! Note that CG modifies the array of the electron mass. 63! 64! Riccardo Bertossa, Federico Grasselli 65! (SISSA - via Bonomea, 265 - 34136 Trieste ITALY) 66!******************************************************************************** 67!06/2018 68!Implemented on the fly change of timestep (with correct rescaling of the 69!velocities) 70! 71! Riccardo Bertossa 72!******************************************************************************** 73 74 75 76MODULE cp_autopilot 77 !--------------------------------------------------------------------------- 78 ! 79 ! This module handles the Autopilot Feature Suite 80 ! Written by Lee Atkinson, with help from the ATP team at Targacept, Inc 81 ! Created June 2005 82 ! Modified by Yonas Abrahm Sept 2006 83 ! 84 ! The address for Targacept, Inc. is: 85 ! 200 East First Street, Suite 86 ! 300, Winston-Salem, North Carolina 27101; 87 ! Attn: Molecular Design. 88 ! 89 ! See README.AUTOPILOT in the Doc directory for more information. 90 !--------------------------------------------------------------------------- 91 92 USE kinds 93 USE autopilot, ONLY : current_nfi, pilot_p, pilot_unit, pause_p,auto_error, & 94 & parse_mailbox, rule_isave, rule_iprint, rule_dt, rule_emass, & 95 & rule_electron_dynamics, rule_electron_damping, rule_ion_dynamics, & 96 & rule_ion_damping, rule_ion_temperature, rule_tempw, & 97 & rule_electron_orthogonalization, rule_tprint, & 98 & rule_nhpcl, rule_fnosep 99 USE autopilot, ONLY : event_index, event_step, event_isave, event_iprint, & 100 & event_dt, event_emass, event_electron_dynamics, event_electron_damping, & 101 & event_ion_dynamics, event_ion_damping, event_ion_temperature, event_tempw, & 102 & event_electron_orthogonalization, & 103 & event_tprint, & 104 & event_nhpcl, event_fnosep 105 106 IMPLICIT NONE 107 SAVE 108 109 110 PRIVATE 111 PUBLIC :: pilot, employ_rules 112 LOGICAL, PRIVATE :: had_tcg_true = .false. 113 LOGICAL, PRIVATE :: had_tens_true = .false. 114CONTAINS 115 116 117 118 119 !----------------------------------------------------------------------- 120 ! EMPLOY_RULES 121 !----------------------------------------------------------------------- 122 SUBROUTINE employ_rules() 123 USE input_parameters, ONLY : dt, & 124 & electron_dynamics, electron_damping, & 125 & ion_dynamics, ion_damping, & 126 & ion_temperature, fnosep, nhpcl, nhptyp, nhgrp, fnhscl, ndega, nat, & 127 & orthogonalization 128 use ions_nose, ONLY: tempw 129 USE control_flags, only: tsde, tsdp, tfor, tcp, tnosep, isave,iprint,& 130 tconvthrs, tolp, & 131 ekin_conv_thr, forc_conv_thr, etot_conv_thr,& 132 tortho, tfirst, tlast, tprint 133 use wave_base, only: frice 134 use ions_base, only: fricp 135 USE ions_nose, ONLY: ions_nose_init,xnhp0, xnhpm, ions_nose_deallocate 136 USE ions_positions, ONLY : taus, tausm 137 USE io_global, ONLY: ionode, ionode_id 138 USE time_step, ONLY : set_time_step,tps 139 USE cp_electronic_mass, ONLY: emass, emass_cutoff, emass_precond 140 USE cg_module, ONLY : tcg,allocate_cg,cg_info, & 141 nfi_firstcg,c0old 142 USE wavefunctions, ONLY : cm_bgrp,c0_bgrp 143 USE ensemble_dft, ONLY : tens,allocate_ensemble_dft 144 USE uspp, ONLY : nkb, nkbus 145 USE electrons_base, ONLY : nspin, nbsp, nbspx, nudx 146 USE gvecw, ONLY : ngw 147 USE fft_base, ONLY : dffts 148 USE cp_main_variables, ONLY : idesc, lambdap, lambda, lambdam 149 USE ions_base, ONLY : nat_ions_base => nat 150 USE cell_base, ONLY : tpiba2 151 USE cp_main_variables, ONLY : ema0bg 152 USE gvecw, ONLY : g2kin, ngw 153 154 IMPLICIT NONE 155 156 REAL(DP) :: old_tps 157 158 !---------------------------------------- 159 ! &CONTROL 160 !---------------------------------------- 161 162 ! ISAVE 163 if (event_isave(event_index)) then 164 isave = rule_isave(event_index) 165 IF ( ionode ) write(*,'(4X,A,15X,I10)') 'Rule event: isave', isave 166 endif 167 168 ! IPRINT 169 if (event_iprint(event_index)) then 170 iprint = rule_iprint(event_index) 171 IF ( ionode ) write(*,'(4X,A,13X,I10)') 'Rule event: iprint', iprint 172 endif 173 174 ! TPRINT 175 if (event_tprint(event_index)) then 176 tprint = rule_tprint(event_index) 177 IF ( ionode ) write(*,*) 'NUOVO Rule event: tprint', tprint 178 endif 179 180 ! DT 181 if (event_dt(event_index)) then 182 IF ( ionode ) write(*,'(4X,A,18X,F10.4)') 'Rule event: dt',rule_dt(event_index) 183 IF (dt /= rule_dt(event_index) ) THEN 184 IF (.NOT. ( tcg .OR. tsde ) ) THEN !if I am doing verlet for electrons 185 !rescale electron velocities (wavefunctions and lagrange 186 !multipliers) 187 IF ( ionode ) write(*,*) 'Rescaling electronic velocities with new dt' 188 lambdam = lambda - (lambda-lambdam)*rule_dt(event_index)/dt 189 cm_bgrp = c0_bgrp - (c0_bgrp-cm_bgrp)*rule_dt(event_index)/dt 190 END IF 191 IF ( (.NOT. tsdp) .AND. tfor ) THEN !if I am doing verlet for ions 192 !rescale ionic velocities (ions and nose variables) 193 IF ( ionode ) write(*,*) 'Rescaling ionic velocities with new dt' 194 tausm = taus - (taus-tausm)*rule_dt(event_index)/dt 195 xnhpm = xnhp0 - (xnhp0-xnhpm)*rule_dt(event_index)/dt 196 END IF 197 END IF 198 dt = rule_dt(event_index) 199 old_tps=tps 200 CALL set_time_step( dt ) !this subroutine sets to zero the elapsed time tps 201 tps=old_tps 202 endif 203 204 !---------------------------------------- 205 ! &SYSTEM 206 !---------------------------------------- 207 208 !---------------------------------------- 209 ! &ELECTRONS 210 !---------------------------------------- 211 212 213 if (event_electron_orthogonalization(event_index)) then 214 orthogonalization=rule_electron_orthogonalization(event_index) 215 select case (orthogonalization) 216 case ('ORTHO') 217 tortho=.true. 218 IF ( ionode ) write(*,*) 'Wow, setting tortho=.true. !' 219 case ('GRAM-SCHMIDT') 220 tortho=.false. 221 IF ( ionode ) write(*,*) 'Wow, setting tortho=.false. ! (Ma cossa xe sta monada?)' 222 case default 223 call auto_error(' autopilot ',' unknown orthogonalization'//trim(orthogonalization) ) 224 end select 225 226 endif 227 228 ! EMASS 229 if (event_emass(event_index)) then 230 emass = rule_emass(event_index) 231 IF ( ionode ) write(*,'(4X,A,15X,F10.4)') 'Rule event: emass', emass 232 endif 233 234 ! ELECTRON_DYNAMICS 235 ! electron_dynamics = 'sd' | 'verlet' | 'damp' | 'none' 236 if (event_electron_dynamics(event_index)) then 237 electron_dynamics= rule_electron_dynamics(event_index) 238 frice = 0.d0 239 ! the cg algorithm uses cm for its internal purposes. The true cm is 240 ! c0old, so now I have to copy it into cm 241 if ( tcg ) then 242 cm_bgrp=c0old 243 ! tfirst=.true. !check if this is needed (I don't think so) 244 ! the conjugate gradient method modifies the electron mass mu(k), 245 ! by calling the routine 246 ! emass_precond_tpa( ema0bg, tpiba2, emass_cutoff ) 247 ! so I call here the standard routine to set mu(k) 248 CALL emass_precond( ema0bg, g2kin, ngw, tpiba2, emass_cutoff) 249 250 endif 251 select case ( electron_dynamics ) 252 case ('SD') 253 tsde = .true. 254 tcg=.false. 255 case ('VERLET') 256 tsde = .false. 257 tcg=.false. 258 case ('DAMP') 259 tsde = .false. 260 frice = electron_damping 261 tcg=.false. 262 case ('NONE') 263 tsde = .false. 264 tcg=.false. 265 case ('CG') 266 tsde=.false. 267 IF ( ionode ) write(*,*) 'Wow, setting tcg=.true. at step '& 268 ,current_nfi,' ! (La ghe domandi a mia molie)' 269 if (.not. tcg) then 270 nfi_firstcg=current_nfi !the first step is different! 271 end if 272 if (.not. had_tcg_true) then 273 had_tcg_true=.true. 274 if (.not. had_tens_true) then 275 CALL allocate_ensemble_dft( nkb, nbsp, ngw, nudx, nspin, nbspx, & 276 dffts%nnr, nat_ions_base, idesc ) 277 endif 278 CALL allocate_cg( ngw, nbspx,nkbus ) 279 endif 280 CALL cg_info() 281 tcg = .true. 282 case default 283 call auto_error(' autopilot ',' unknown electron_dynamics '//trim(electron_dynamics) ) 284 end select 285 286 IF ( ionode ) write(*,'(4X,A,2X,A10)') 'Rule event: electron_dynamics', electron_dynamics 287 288 endif 289 290 291 ! ELECTRON_DAMPING 292 if (event_electron_damping(event_index)) then 293 ! meaningful only if " electron_dynamics = 'damp' " 294 electron_damping = rule_electron_damping(event_index) 295 frice = electron_damping 296 IF ( ionode ) write(*,'(4X,A,4X,F10.4)') 'Rule event: electron_damping', electron_damping 297 endif 298 299 !---------------------------------------- 300 ! &IONS 301 !---------------------------------------- 302 303 304 ! ION_DYNAMICS 305 ! ion_dynamics = 'sd' | 'verlet' | 'damp' | 'none' 306 if (event_ion_dynamics(event_index)) then 307 ion_dynamics= rule_ion_dynamics(event_index) 308 tconvthrs%active = .FALSE. 309 tconvthrs%nstep = 1 310 tconvthrs%ekin = 0.0d0 311 tconvthrs%derho = 0.0d0 312 tconvthrs%force = 0.0d0 313 314 select case ( ion_dynamics ) 315 case ('SD') 316 tsdp = .true. 317 tfor = .true. 318 fricp= 0.d0 319 tconvthrs%ekin = ekin_conv_thr 320 tconvthrs%derho = etot_conv_thr 321 tconvthrs%force = forc_conv_thr 322 tconvthrs%active = .TRUE. 323 tconvthrs%nstep = 1 324 case ('VERLET') 325 tsdp = .false. 326 tfor = .true. 327 fricp= 0.d0 328 case ('DAMP') 329 tsdp = .false. 330 tfor = .true. 331 fricp= ion_damping 332 tconvthrs%ekin = ekin_conv_thr 333 tconvthrs%derho = etot_conv_thr 334 tconvthrs%force = forc_conv_thr 335 tconvthrs%active = .TRUE. 336 tconvthrs%nstep = 1 337 case ('NONE') 338 tsdp = .false. 339 tfor = .false. 340 fricp= 0.d0 341 case default 342 call auto_error(' iosys ',' unknown ion_dynamics '//trim(ion_dynamics) ) 343 end select 344 345 endif 346 347 348 ! ION_DAMPING 349 if (event_ion_damping(event_index)) then 350 ! meaningful only if " ion_dynamics = 'damp' " 351 ion_damping = rule_ion_damping(event_index) 352 IF ( ionode ) write(*,'(4X,A,9X,F10.4)') 'Rule event: ion_damping', ion_damping 353 endif 354 355 356 ! ION_TEMPERATURE 357 if (event_ion_temperature(event_index)) then 358 ion_temperature = rule_ion_temperature(event_index) 359 tcp = .FALSE. 360 tnosep = .FALSE. 361 tolp = tolp 362 select case ( ion_temperature ) 363 ! temperature control of ions via nose' thermostat 364 ! tempw (real(DP)) frequency (in which units?) 365 ! fnosep (real(DP)) temperature (in which units?) 366 case ('NOSE') 367 tnosep = .true. 368 tcp = .false. 369 if ( .not. event_tempw(event_index)) then 370 IF ( ionode ) write(*,*) 'WARNING: missing tempw event (if undefined can make the code bananas) ' 371 end if 372 if ( .not. tfor) then 373 IF ( ionode ) write(*,*) 'WARNING: not doing Verlet on ions? (am I correct?)' 374 end if 375 case ('NOT_CONTROLLED') 376 tnosep = .false. 377 tcp = .false. 378 case ('RESCALING' ) 379 tnosep = .false. 380 tcp = .true. 381 case default 382 call auto_error(' iosys ',' unknown ion_temperature '//trim(ion_temperature) ) 383 end select 384 385 IF ( ionode ) write(*,'(4X,A,5X,A)') 'Rule event: ion_temperature', ion_temperature 386 387 endif 388 389 ! TEMPW 390 if (event_tempw(event_index)) then 391 tempw = rule_tempw(event_index) 392 ! The follwiong is a required side effect 393 ! when resetting tempw 394 IF ( ionode ) write(*,'(4X,A,15X,F10.4)') 'Rule event: tempw', tempw 395 endif 396 397 ! NHPCL 398 if (event_nhpcl(event_index)) then 399 nhpcl = rule_nhpcl(event_index) 400 call ions_nose_deallocate() ! the dimension of the nose array can change 401 if ( ionode ) write(*,'(4X,A,15X,I1)') 'Rule event: nhpcl', nhpcl 402 403 endif 404 405 ! FNOSEP 406 ! TODO: implement reading of full FNOSEP array 407 if (event_fnosep(event_index)) then 408 fnosep = rule_fnosep(event_index) 409 if ( ionode ) then 410 write(*,'(4X,A,15X,F10.4)') 'Rule event: fnosep', fnosep(1) 411 write(*,'(4X,A)') 'WARNING: all fnosep of the chain are set to the same value' 412 endif 413 414 !IF ( ionode ) write(*,*) 'Rule event: fnosep', fnosep 415 endif 416 417 ! nose initialization 418 if (event_tempw(event_index) .or. event_nhpcl(event_index) .or. event_fnosep(event_index) .or. & 419 (event_ion_temperature(event_index) .and. rule_ion_temperature(event_index) .eq. 'NOSE' ) & 420 ) then 421 call ions_nose_init( tempw, fnosep, nhpcl, nhptyp, ndega, nhgrp, fnhscl) 422 endif 423 424 !---------------------------------------- 425 ! &CELL 426 !---------------------------------------- 427 428 !---------------------------------------- 429 ! &PHONON 430 !---------------------------------------- 431 432 433 END SUBROUTINE employ_rules 434 435 436 437 !----------------------------------------------------------------------- 438 ! PILOT 439 ! 440 ! Here is the main pilot routine called in CPR, at the top 441 ! of the basic dynamics loop just after nose hoover update 442 !----------------------------------------------------------------------- 443 subroutine pilot (nfi) 444 USE parser, ONLY: parse_unit 445 USE io_global, ONLY: ionode, ionode_id 446 USE mp, ONLY : mp_bcast, mp_barrier 447 USE mp_world, ONLY : world_comm 448#if defined (__NAG) 449 USE f90_unix_proc 450#endif 451 USE ensemble_dft, ONLY : tens 452 USE cg_module, ONLY : tcg 453 USE control_flags , only : tprint 454 455 IMPLICIT NONE 456 INTEGER :: nfi 457 LOGICAL :: file_p 458 CHARACTER (LEN=256) :: mbfile = "pilot.mb" 459 460 ! The tcg control variable, that enables conjugate gradient electron 461 ! dynamics, allocates some arrays when used in the input file. So I check 462 ! if this variable is defined, and set had_tcg_true to true. 463 464 if (tcg) had_tcg_true = .true. 465 if (tens) had_tens_true = .true. 466 467 ! Dynamics Loop Started 468 pilot_p = .TRUE. 469 470 ! This is so we can usurp the exiting parser 471 ! that defaults to stdin (unit=5) 472 ! We have to do it this way if we are to 473 ! call (reuse) the card_autopilot that is called 474 ! by read_cards 475 parse_unit = pilot_unit 476 477 ! Our own local for nfi 478 current_nfi = nfi 479 480 481 ! This allows one pass. Calling parse_mailbox will either: 482 ! 1) call init_auto_pilot, which will always set this modules global PAUSE_P variable to FALSE 483 ! 2) detect a pause indicator, setting PAUSE_P to TRUE until a new mailbox overrides. 484 pause_loop: do 485 486 file_p = .FALSE. 487 IF ( ionode ) INQUIRE( FILE = TRIM( mbfile ), EXIST = file_p ) 488 call mp_bcast(file_p, ionode_id,world_comm) 489 490 IF ( file_p ) THEN 491 492 IF ( ionode ) THEN 493 WRITE(*,*) 494 WRITE(*,*) '****************************************************' 495 WRITE(*,*) ' Autopilot: Mailbox found at nfi=', current_nfi 496 END IF 497 FLUSH(6) 498 499 ! Open the mailbox 500 IF ( ionode ) OPEN( UNIT = pilot_unit, FILE = TRIM( mbfile ) ) 501 502 ! Will reset PAUSE_P to false unless there is a PAUSE cmd 503 ! The following call is MPI safe! It only generates side effects 504 CALL parse_mailbox() 505 call mp_barrier( world_comm ) 506 507 IF ( ionode ) THEN 508 WRITE(*,*) ' Autopilot: Done reading mailbox' 509 WRITE(*,*) '****************************************************' 510 WRITE(*,*) 511 END IF 512 513 ! Perhaps instead of deleting move the file as an input log 514 IF( ionode ) CLOSE( UNIT = pilot_unit, STATUS = 'DELETE' ) 515 516 END IF 517 518 IF( .NOT. pause_p ) THEN 519 EXIT pause_loop 520 ELSE 521 IF( ionode ) write(*,*) 'SLEEPING .... send another pilot.mb' 522 call sleep (5) 523 END if 524 525 end do pause_loop 526 527 528 529 530 ! Autopilot (Dynamic Rules) Implementation 531 ! When nfi has passed (is greater than 532 ! the next event, then employ rules 533 ! Mailbox may have issued several rules 534 ! Attempt to catch up! 535 do while (current_nfi >= event_step(event_index) ) 536 537 IF ( ionode ) THEN 538 WRITE(*,*) 539 WRITE(*,*) '****************************************************' 540 WRITE(*,*) ' Autopilot employ rules: ' 541 END IF 542 call employ_rules() 543 IF ( ionode ) THEN 544 WRITE(*,*) '****************************************************' 545 WRITE(*,*) 546 END IF 547 call mp_barrier( world_comm ) 548 549 ! update event_index to current 550 event_index = event_index + 1 551 552 enddo 553 554 ! During the last step of electron conjugate gradient, 555 ! set tprint to .true. if next step has electron_dynamics = Verlet 556 ! This is needed to calculate wavefunctions at time t and (t - dt) 557 if (need_tprint_true() ) then 558 tprint = .true. 559 if (ionode) then 560 WRITE(*,*) '==================================================' 561 WRITE(*,*) ' Setting tprint=.true. for this step (last of CG)' 562 WRITE(*,*) '==================================================' 563 endif 564 endif 565 566 end subroutine pilot 567 568 function need_tprint_true() 569 ! Check if next step is a Verlet. In such case returns .true. 570 ! Used whenever tprint == .true. is needed, e.g. to let CG 571 ! calculate wavefunctions at time t and (t - dt) via projections 572 ! onto the occupied manifold. 573 ! This function has to be called after 'call employ_rules()'! 574 USE cg_module, ONLY : tcg 575 LOGICAL :: need_tprint_true 576 INTEGER :: event_idx 577 578 need_tprint_true = .FALSE. 579 event_idx = event_index 580 do while ( event_idx .le. size(event_step) .and. (event_step(event_idx)==current_nfi+1) ) 581 if ( tcg .and. event_electron_dynamics(event_idx) .and. & 582 & (rule_electron_dynamics(event_idx)=='VERLET') ) then 583 need_tprint_true = .TRUE. 584 end if 585 event_idx = event_idx + 1 586 enddo 587 RETURN 588 end function need_tprint_true 589 590END MODULE cp_autopilot 591 592