1! 2! Copyright (C) 2001-2011 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! TB 10! included test if relaxz=.TRUE. to allow for movement of the center of mass 11! search for 'TB' 12!---------------------------------------------------------------------------- 13! 14#undef __NPT 15#if defined (__NPT) 16#define RELAXTIME 2000.D0 17#define TARGPRESS 2.39D0 18#endif 19! 20!---------------------------------------------------------------------------- 21MODULE dynamics_module 22 !---------------------------------------------------------------------------- 23 !! It includes all the quantities and procedures related to the molecular 24 !! dynamics. 25 ! 26 USE kinds, ONLY : DP 27 USE ions_base, ONLY : amass 28 USE io_global, ONLY : stdout 29 USE io_files, ONLY : prefix, tmp_dir, seqopn 30 USE constants, ONLY : tpi, fpi 31 USE constants, ONLY : amu_ry, ry_to_kelvin, au_ps, bohr_radius_cm, & 32 ry_kbar 33 USE constants, ONLY : eps8 34 USE control_flags, ONLY : tolp 35 ! 36 USE basic_algebra_routines 37 ! 38 IMPLICIT NONE 39 ! 40 SAVE 41 PRIVATE 42 PUBLIC :: verlet, proj_verlet, terminate_verlet, & 43 langevin_md, smart_MC, allocate_dyn_vars, deallocate_dyn_vars 44 PUBLIC :: temperature, refold_pos, vel 45 PUBLIC :: dt, delta_t, nraise, control_temp, thermostat 46 ! 47 ! 48 REAL(DP) :: dt 49 !! time step 50 REAL(DP) :: temperature 51 !! starting temperature 52 REAL(DP) :: virial 53 !! virial (used for the pressure) 54 REAL(DP) :: delta_t 55 !! parameter used in thermalization 56 INTEGER :: nraise 57 !! parameter used in thermalization 58 INTEGER :: ndof 59 !! the number of degrees of freedom 60 INTEGER :: num_accept=0 61 !! Number of the accepted proposal in Smart_MC 62 LOGICAL :: vel_defined 63 !! if true, vel is used rather than tau_old to do the next step 64 LOGICAL :: control_temp 65 !! if true a thermostat is used to control the temperature 66 LOGICAL :: refold_pos 67 !! if true the positions are refolded into the supercell 68 LOGICAL :: first_iter=.TRUE. 69 !! true if this is the first ionic iteration 70 CHARACTER(len=10) :: thermostat 71 !! the thermostat used to control the temperature 72 REAL(DP), ALLOCATABLE :: tau_smart(:,:) 73 !! used for smart Monte Carlo to store the atomic position of the 74 !! previous step. 75 REAL(DP), ALLOCATABLE :: force_smart(:,:) 76 !! used for smart Monte Carlo to store the force of the 77 !! previous step. 78 REAL(DP) :: etot_smart 79 !! used to keep the energy of the previous iteration 80 REAL(DP), ALLOCATABLE :: tau_old(:,:) 81 !! the atomic positions at the previous iteration 82 REAL(DP), ALLOCATABLE :: tau_new(:,:) 83 !! the atomic positions at the new iteration 84 REAL(DP), ALLOCATABLE :: tau_ref(:,:) 85 !! reference atomic positions 86 REAL(DP), ALLOCATABLE :: vel(:,:) 87 !! velocities 88 REAL(DP), ALLOCATABLE :: acc(:,:) 89 !! accelerations 90 REAL(DP), ALLOCATABLE :: chi(:,:) 91 !! chi 92 REAL(DP), ALLOCATABLE :: mass(:) 93 !! atomic masses 94 REAL(DP), ALLOCATABLE :: diff_coeff(:) 95 !! diffusion coefficients 96 REAL(DP), ALLOCATABLE :: radial_distr(:,:) 97 !! radial distribution 98 ! 99 INTEGER, PARAMETER :: hist_len = 1000 100 ! 101CONTAINS 102 ! 103 ! ... public methods 104 ! 105 !------------------------------------------------------------------------ 106 SUBROUTINE allocate_dyn_vars() 107 !------------------------------------------------------------------------ 108 !! Allocates dynamics variables 109 ! 110 USE ions_base, ONLY : nat 111 ! 112 IF ( .NOT.ALLOCATED( mass ) ) ALLOCATE( mass( nat ) ) 113 ! 114 IF ( .NOT.ALLOCATED( tau_old ) ) ALLOCATE( tau_old( 3, nat ) ) 115 IF ( .NOT.ALLOCATED( tau_new ) ) ALLOCATE( tau_new( 3, nat ) ) 116 IF ( .NOT.ALLOCATED( tau_ref ) ) ALLOCATE( tau_ref( 3, nat ) ) 117 ! 118 IF ( .NOT.ALLOCATED( vel ) ) ALLOCATE( vel( 3, nat ) ) 119 IF ( .NOT.ALLOCATED( acc ) ) ALLOCATE( acc( 3, nat ) ) 120 IF ( .NOT.ALLOCATED( chi ) ) ALLOCATE( chi( 3, nat ) ) 121 ! 122 IF ( .NOT.ALLOCATED( diff_coeff ) ) ALLOCATE( diff_coeff( nat ) ) 123 ! 124 IF ( .NOT.ALLOCATED( radial_distr ) ) & 125 ALLOCATE( radial_distr( hist_len , nat ) ) 126 ! 127 END SUBROUTINE allocate_dyn_vars 128 ! 129 ! 130 !------------------------------------------------------------------------ 131 SUBROUTINE deallocate_dyn_vars() 132 !------------------------------------------------------------------------ 133 !! Deallocates dynamics variables. 134 ! 135 IF ( ALLOCATED( mass ) ) DEALLOCATE( mass ) 136 IF ( ALLOCATED( tau_old ) ) DEALLOCATE( tau_old ) 137 IF ( ALLOCATED( tau_new ) ) DEALLOCATE( tau_new ) 138 IF ( ALLOCATED( tau_ref ) ) DEALLOCATE( tau_ref ) 139 IF ( ALLOCATED( vel ) ) DEALLOCATE( vel ) 140 IF ( ALLOCATED( acc ) ) DEALLOCATE( acc ) 141 IF ( ALLOCATED( chi ) ) DEALLOCATE( chi ) 142 IF ( ALLOCATED( diff_coeff ) ) DEALLOCATE( diff_coeff ) 143 IF ( ALLOCATED( radial_distr ) ) DEALLOCATE( radial_distr ) 144 ! 145 END SUBROUTINE deallocate_dyn_vars 146 ! 147 ! 148 !------------------------------------------------------------------------ 149 SUBROUTINE verlet() 150 !------------------------------------------------------------------------ 151 !! This routine performs one step of molecular dynamics evolution 152 !! using the Verlet algorithm. 153 !! The parameters: 154 ! 155 !! * mass: mass of the atoms; 156 !! * dt: time step; 157 !! * temperature: starting temperature. 158 ! 159 !! The starting velocities of atoms are set accordingly 160 !! to the starting temperature, in random directions. 161 !! The initial velocity distribution is therefore a 162 !! constant. 163 ! 164 !! Dario Alfe' 1997 and Carlo Sbraccia 2004-2006. 165 ! 166 USE ions_base, ONLY : nat, nsp, ityp, tau, if_pos, atm 167 USE cell_base, ONLY : alat, omega 168 USE ener, ONLY : etot 169 USE force_mod, ONLY : force, lstres 170 USE control_flags, ONLY : istep, lconstrain, tv0rd 171 ! 172 USE constraints_module, ONLY : nconstr, check_constraint 173 USE constraints_module, ONLY : remove_constr_force, remove_constr_vec 174 ! 175 IMPLICIT NONE 176 ! 177 ! ... local variables 178 ! 179 REAL(DP) :: ekin, etotold 180 REAL(DP) :: total_mass, temp_new, temp_av, elapsed_time 181 REAL(DP) :: delta(3), ml(3), mlt 182 INTEGER :: na 183 ! istep counts all MD steps, including those of previous runs 184#if defined (__NPT) 185 REAL(DP) :: chi, press_new 186#endif 187 LOGICAL :: file_exists, leof 188 REAL(DP), EXTERNAL :: dnrm2 189 REAL(DP) :: kstress(3,3) 190 INTEGER :: i, j 191 ! 192 ! ... the number of degrees of freedom 193 ! 194 IF ( ANY( if_pos(:,:) == 0 ) ) THEN 195 ! 196 ndof = 3*nat - COUNT( if_pos(:,:) == 0 ) - nconstr 197 ! 198 ELSE 199 ! 200 ndof = 3*nat - 3 - nconstr 201 ! 202 ENDIF 203 ! 204 vel_defined = .TRUE. 205 temp_av = 0.D0 206 ! 207 CALL seqopn( 4, 'md', 'FORMATTED', file_exists ) 208 ! 209 IF ( file_exists ) THEN 210 ! 211 ! ... the file is read : simulation is continuing 212 ! 213 READ( UNIT = 4, FMT = * ) etotold, istep, tau_old(:,:), leof 214 ! 215 IF ( leof ) THEN 216 ! 217 ! ... the file was created by projected_verlet: Ignore it 218 ! 219 CALL md_init() 220 ! 221 ELSE 222 ! 223 vel_defined = .FALSE. 224 ! 225 READ( UNIT = 4, FMT = * ) & 226 temp_new, temp_av, mass(:), total_mass, elapsed_time, & 227 tau_ref(:,:) 228 ! 229 ENDIF 230 ! 231 CLOSE( UNIT = 4, STATUS = 'KEEP' ) 232 ! 233 ELSE 234 ! 235 CLOSE( UNIT = 4, STATUS = 'DELETE' ) 236 ! 237 ! ... the file is absent : simulation is starting from scratch 238 ! 239 CALL md_init() 240 ! 241 ENDIF 242 ! 243 ! ... elapsed_time is in picoseconds 244 ! 245 elapsed_time = elapsed_time + dt*2.D0*au_ps 246 ! 247 istep = istep + 1 248 ! 249 WRITE( UNIT = stdout, & 250 FMT = '(/,5X,"Entering Dynamics:",T28,"iteration",T37," = ", & 251 &I5,/,T28,"time",T37," = ",F8.4," pico-seconds",/)' ) & 252 istep, elapsed_time 253 ! 254 IF ( control_temp ) CALL apply_thermostat() 255 ! 256 ! ... we first remove the component of the force along the 257 ! ... constraint gradient ( this constitutes the initial 258 ! ... guess for the calculation of the lagrange multipliers ) 259 ! 260 IF ( lconstrain ) & 261 CALL remove_constr_force( nat, tau, if_pos, ityp, alat, force ) 262 ! 263 ! ... calculate accelerations in a.u. units / alat 264 ! 265 FORALL( na = 1:nat ) acc(:,na) = force(:,na) / mass(na) / alat 266 ! 267 ! ... Verlet integration scheme 268 ! 269 IF (vel_defined) THEN 270 ! 271 ! ... remove the component of the velocity along the 272 ! ... constraint gradient 273 ! 274 IF ( lconstrain ) & 275 CALL remove_constr_vec( nat, tau, if_pos, ityp, alat, vel ) 276 ! 277 tau_new(:,:) = tau(:,:) + vel(:,:) * dt + 0.5_DP * acc(:,:) * dt**2 278 tau_old(:,:) = tau(:,:) - vel(:,:) * dt + 0.5_DP * acc(:,:) * dt**2 279 ! 280 ELSE 281 ! 282 tau_new(:,:) = 2.D0*tau(:,:) - tau_old(:,:) + acc(:,:) * dt**2 283 ! 284 ENDIF 285 ! 286 IF ( .NOT. ANY( if_pos(:,:) == 0 ) ) THEN 287 ! 288 ! ... if no atom has been fixed we compute the displacement of the 289 ! ... center of mass and we subtract it from the displaced positions 290 ! 291 delta(:) = 0.D0 292 DO na = 1, nat 293 delta(:) = delta(:) + mass(na)*( tau_new(:,na) - tau(:,na) ) 294 ENDDO 295 delta(:) = delta(:) / total_mass 296 FORALL( na = 1:nat ) tau_new(:,na) = tau_new(:,na) - delta(:) 297 ! 298 IF (vel_defined) THEN 299 delta(:) = 0.D0 300 DO na = 1, nat 301 delta(:) = delta(:) + mass(na)*( tau_old(:,na) - tau(:,na) ) 302 ENDDO 303 delta(:) = delta(:) / total_mass 304 FORALL( na = 1:nat ) tau_old(:,na) = tau_old(:,na) - delta(:) 305 ENDIF 306 ! 307 ENDIF 308 ! 309 IF ( lconstrain ) THEN 310 ! 311 ! ... check if the new positions satisfy the constrain equation 312 ! 313 CALL check_constraint( nat, tau_new, tau, & 314 force, if_pos, ityp, alat, dt, amu_ry ) 315 ! 316#if ! defined (__REDUCE_OUTPUT) 317 ! 318 WRITE( stdout, '(/,5X,"Constrained forces (Ry/au):",/)') 319 ! 320 DO na = 1, nat 321 ! 322 WRITE( stdout, & 323 '(5X,"atom ",I3," type ",I2,3X,"force = ",3F14.8)' ) & 324 na, ityp(na), force(:,na) 325 ! 326 ENDDO 327 ! 328 WRITE( stdout, '(/5X,"Total force = ",F12.6)') dnrm2( 3*nat, force, 1 ) 329 ! 330#endif 331 332 IF (vel_defined) THEN 333 CALL check_constraint( nat, tau_old, tau, & 334 force, if_pos, ityp, alat, dt, amu_ry ) 335 ENDIF 336 ! 337 ENDIF 338 ! 339 ! ... the linear momentum and the kinetic energy are computed here 340 ! 341 vel = ( tau_new - tau_old ) / ( 2.D0*dt ) * DBLE( if_pos ) 342 ! 343 ml = 0.D0 344 ekin = 0.D0 345 kstress = 0.d0 346 ! 347 DO na = 1, nat 348 ! 349 ml(:) = ml(:) + vel(:,na) * mass(na) 350 ekin = ekin + 0.5D0 * mass(na) * & 351 ( vel(1,na)**2 + vel(2,na)**2 + vel(3,na)**2 ) 352 DO i = 1, 3 353 DO j = 1, 3 354 kstress(i,j) = kstress(i,j) + mass(na)*vel(i,na)*vel(j,na) 355 ENDDO 356 ENDDO 357 ! 358 ENDDO 359 ! 360 ekin = ekin*alat**2 361 kstress = kstress * alat**2 / omega 362 ! 363 ! ... find the new temperature and update the average 364 ! 365 temp_new = 2.D0 / DBLE( ndof ) * ekin * ry_to_kelvin 366 ! 367 temp_av = temp_av + temp_new 368 ! 369#if defined (__NPT) 370 ! 371 ! ... find the new pressure (in Kbar) 372 ! 373 press_new = ry_kbar*( nat*temp_new/ry_to_kelvin + virial ) / omega 374 ! 375 chi = 1.D0 - dt / RELAXTIME*( TARGPRESS - press_new ) 376 ! 377 omega = chi * omega 378 alat = chi**(1.D0/3.D0) * alat 379 ! 380 WRITE( stdout, '(/,5X,"NEW ALAT = ",F8.5,2X,"Bohr" )' ) alat 381 WRITE( stdout, '( 5X,"PRESSURE = ",F8.5,2X,"Kbar",/)' ) press_new 382 ! 383#endif 384 ! 385 ! ... save all the needed quantities on file 386 ! 387 CALL seqopn( 4, 'md', 'FORMATTED', file_exists ) 388 ! 389 leof = .FALSE. 390 WRITE( UNIT = 4, FMT = * ) etot, istep, tau(:,:), leof 391 ! 392 WRITE( UNIT = 4, FMT = * ) & 393 temp_new, temp_av, mass(:), total_mass, elapsed_time, tau_ref(:,:) 394 ! 395 CLOSE( UNIT = 4, STATUS = 'KEEP' ) 396 ! 397 ! ... here the tau are shifted 398 ! 399 tau(:,:) = tau_new(:,:) 400 ! 401#if ! defined (__REDUCE_OUTPUT) 402 ! 403 CALL output_tau( .FALSE., .FALSE. ) 404 ! 405#endif 406 ! 407 ! ... infos are written on the standard output 408 ! 409 WRITE( stdout, '(5X,"kinetic energy (Ekin) = ",F20.8," Ry",/, & 410 & 5X,"temperature = ",F20.8," K ",/, & 411 & 5X,"Ekin + Etot (const) = ",F20.8," Ry")' ) & 412 ekin, temp_new, ( ekin + etot ) 413 ! 414 IF (lstres) WRITE ( stdout, & 415 '(5X,"Ions kinetic stress = ",F15.2," (kbar)",/3(27X,3F15.2/)/)') & 416 ((kstress(1,1)+kstress(2,2)+kstress(3,3))/3.d0*ry_kbar), & 417 (kstress(i,1)*ry_kbar,kstress(i,2)*ry_kbar,kstress(i,3)*ry_kbar, i=1,3) 418 ! 419 IF ( .NOT.( lconstrain .or. ANY( if_pos(:,:) == 0 ) ) ) THEN 420 ! 421 ! ... total linear momentum must be zero if all atoms move 422 ! 423 mlt = norm( ml(:) ) 424 ! 425 IF ( mlt > eps8 ) & 426 CALL infomsg( 'dynamics', 'Total linear momentum <> 0' ) 427 ! 428 WRITE( stdout, '(/,5X,"Linear momentum :",3(2X,F14.10))' ) ml(:) 429 ! 430 ENDIF 431 ! 432 ! ... compute the average quantities 433 ! 434 CALL compute_averages( istep ) 435 ! 436 CONTAINS 437 ! 438 !-------------------------------------------------------------------- 439 SUBROUTINE md_init() 440 !-------------------------------------------------------------------- 441 ! 442 IMPLICIT NONE 443 ! 444 istep = 0 445 ! 446 WRITE( UNIT = stdout, & 447 FMT = '(/,5X,"Molecular Dynamics Calculation")' ) 448 ! 449 ! ... atoms are refold in the central box if required 450 ! 451 IF ( refold_pos ) CALL refold_tau() 452 ! 453 ! ... reference positions 454 ! 455 tau_ref(:,:) = tau(:,:) 456 ! 457 IF ( control_temp ) THEN 458 ! 459 WRITE( stdout, & 460 '(/,5X,"Starting temperature",T27," = ",F8.2," K")' ) & 461 temperature 462 ! 463 SELECT CASE( TRIM( thermostat ) ) 464 ! 465 CASE( 'andersen', 'Andersen' ) 466 ! 467 WRITE( UNIT = stdout, & 468 FMT = '(/,5X,"temperature is controlled by Andersen ", & 469 & "thermostat",/,5x,"Collision frequency =",& 470 & f7.4,"/timestep")' ) 1.0_dp/nraise 471 ! 472 CASE( 'berendsen', 'Berendsen' ) 473 ! 474 WRITE( UNIT = stdout, & 475 FMT = '(/,5X,"temperature is controlled by soft ", & 476 & "(Berendsen) velocity rescaling",/,5x,& 477 & "Characteristic time =",i3,"*timestep")') & 478 nraise 479 ! 480 CASE( 'svr', 'Svr', 'SVR' ) 481 ! 482 WRITE( UNIT = stdout, & 483 FMT = '(/,5X,"temperature is controlled by ", & 484 & "stochastic velocity rescaling",/,5x,& 485 & "Characteristic time =",i3,"*timestep")') & 486 nraise 487 CASE( 'initial', 'Initial' ) 488 ! 489 WRITE( UNIT = stdout, & 490 FMT = '(/,5X,"temperature is set once at start"/)' ) 491 ! 492 CASE DEFAULT 493 ! 494 WRITE( UNIT = stdout, & 495 FMT = '(/,5X,"temperature is controlled by ", & 496 & "velocity rescaling (",A,")"/)' )& 497 TRIM( thermostat ) 498 ! 499 END SELECT 500 ! 501 ENDIF 502 ! 503 DO na = 1, nsp 504 ! 505 WRITE( UNIT = stdout, & 506 FMT = '(5X,"mass ",A2,T27," = ",F8.2)' ) atm(na), amass(na) 507 ! 508 ENDDO 509 ! 510 WRITE( UNIT = stdout, & 511 FMT = '(5X,"Time step",T27," = ",F8.2," a.u.,",F8.4, & 512 & " femto-seconds")' ) dt, dt*2.D+3*au_ps 513 ! 514 ! ... masses in rydberg atomic units 515 ! 516 total_mass = 0.D0 517 ! 518 DO na = 1, nat 519 ! 520 mass(na) = amass( ityp(na) ) * amu_ry 521 ! 522 total_mass = total_mass + mass(na) 523 ! 524 ENDDO 525 ! 526 IF ( tv0rd ) THEN ! initial velocities available from input file 527 ! 528 vel(:,:) = vel(:,:) / alat 529 ! 530 ELSEIF ( control_temp ) THEN 531 ! 532 ! ... initial thermalization. N.B. tau is in units of alat 533 ! 534 CALL start_therm() 535 vel_defined = .TRUE. 536 ! 537 temp_new = temperature 538 ! 539 temp_av = 0.D0 540 ! 541 ELSE 542 ! 543 vel(:,:) = 0.0_DP 544 vel_defined = .TRUE. 545 ! 546 ENDIF 547 ! 548 elapsed_time = 0.D0 549 ! 550 END SUBROUTINE md_init 551 ! 552 !-------------------------------------------------------------------- 553 SUBROUTINE apply_thermostat() 554 !-------------------------------------------------------------------- 555 ! 556 USE random_numbers, ONLY : randy, gauss_dist 557 ! 558 IMPLICIT NONE 559 ! 560 INTEGER :: nat_moved 561 REAL(DP) :: sigma, kt 562 ! 563 IF (.NOT. vel_defined) THEN 564 vel(:,:) = (tau(:,:) - tau_old(:,:)) / dt 565 ENDIF 566 ! 567 SELECT CASE( TRIM( thermostat ) ) 568 CASE( 'rescaling' ) 569 ! 570 IF ( ABS(temp_new-temperature) > tolp ) THEN 571 WRITE( UNIT = stdout, & 572 FMT = '(/,5X,"Velocity rescaling: T (",F6.1,"K) ", & 573 & "out of range, reset to " ,F6.1)' ) & 574 temp_new, temperature 575 CALL thermalize( 0, temp_new, temperature ) 576 ENDIF 577 ! 578 CASE( 'rescale-v', 'rescale-V', 'rescale_v', 'rescale_V' ) 579 ! 580 IF ( MOD( istep, nraise ) == 0 ) THEN 581 ! 582 temp_av = temp_av / DBLE( nraise ) 583 ! 584 WRITE( UNIT = stdout, & 585 FMT = '(/,5X,"Velocity rescaling: average T on ",i3, & 586 &" steps (",F6.1,"K) reset to ",F6.1)' ) & 587 nraise, temp_av, temperature 588 ! 589 CALL thermalize( 0, temp_new, temperature ) 590 temp_av = 0.D0 591 ENDIF 592 ! 593 CASE( 'rescale-T', 'rescale-t', 'rescale_T', 'rescale_t' ) 594 ! Clearly it makes sense to check for positive delta_t 595 ! If a negative delta_t is given, I suggest to have a message 596 ! printed, that delta_t is ignored (TODO) 597 IF ( delta_t > 0 ) THEN 598 ! 599 temperature = temperature*delta_t 600 ! 601 WRITE( UNIT = stdout, & 602 FMT = '(/,5X,"Thermalization: T (",F6.1,"K) rescaled ",& 603 & "by a factor ",F6.3)' ) temp_new, delta_t 604 ! 605 CALL thermalize( 0, temp_new, temperature ) 606 ! 607 ENDIF 608 CASE( 'reduce-T', 'reduce-t', 'reduce_T', 'reduce_t' ) 609 IF ( mod( istep, nraise ) == 0 ) THEN 610 ! 611 ! First printing message, than reduce target temperature: 612 ! 613 IF ( delta_t > 0 ) THEN 614 WRITE( UNIT = stdout, & 615 FMT = '(/,5X,"Thermalization: T (",F6.1,"K) augmented ",& 616 & "by ",F6.3)' ) temperature, delta_t 617 618 ELSE 619 WRITE( UNIT = stdout, & 620 FMT = '(/,5X,"Thermalization: T (",F6.1,"K) reduced ",& 621 & "by ",F6.3)' ) temperature, -delta_t 622 ENDIF 623 ! I check whether the temperature is negative, so that I avoid 624 ! nonsensical behavior: 625 IF (temperature < 0.0D0 ) CALL errore('apply_thermostat','Negative target temperature',1) 626 ! 627 temperature = temperature + delta_t 628 ! 629 CALL thermalize( 0, temp_new, temperature ) 630 ! 631 ENDIF 632 ! 633 CASE( 'berendsen', 'Berendsen' ) 634 ! 635 WRITE( UNIT = stdout, & 636 FMT = '(/,5X,"Soft (Berendsen) velocity rescaling")' ) 637 ! 638 CALL thermalize( nraise, temp_new, temperature ) 639 ! 640 CASE( 'svr', 'Svr', 'SVR' ) 641 ! 642 WRITE( UNIT = stdout, & 643 FMT = '(/,5X,"Canonical sampling velocity rescaling")' ) 644 ! 645 CALL thermalize_resamp_vscaling( nraise, temp_new, temperature ) 646 ! 647 CASE( 'andersen', 'Andersen' ) 648 ! 649 kt = temperature / ry_to_kelvin 650 nat_moved = 0 651 ! 652 DO na = 1, nat 653 ! 654 IF ( randy() < 1.D0 / DBLE( nraise ) ) THEN 655 ! 656 nat_moved = nat_moved + 1 657 sigma = SQRT( kt / mass(na) ) 658 ! 659 ! ... N.B. velocities must in a.u. units of alat and are zero 660 ! ... for fixed ions 661 ! 662 vel(:,na) = DBLE( if_pos(:,na) ) * & 663 gauss_dist( 0.D0, sigma, 3 ) / alat 664 ! 665 ENDIF 666 ! 667 ENDDO 668 ! 669 IF ( nat_moved > 0) WRITE( UNIT = stdout, & 670 FMT = '(/,5X,"Andersen thermostat: ",I4," collisions")' ) & 671 nat_moved 672 ! 673 CASE( 'initial', 'Initial' ) 674 ! 675 CONTINUE 676 ! 677 END SELECT 678 ! 679 ! ... the old positions are updated to reflect the new velocities 680 ! 681 IF (.NOT. vel_defined) THEN 682 tau_old(:,:) = tau(:,:) - vel(:,:) * dt 683 ENDIF 684 ! 685 END SUBROUTINE apply_thermostat 686 ! 687 !----------------------------------------------------------------------- 688 SUBROUTINE start_therm() 689 !----------------------------------------------------------------------- 690 !! Starting thermalization of the system. 691 ! 692 USE symm_base, ONLY : invsym, nsym, irt 693 USE cell_base, ONLY : alat 694 USE ions_base, ONLY : nat, if_pos 695 USE random_numbers, ONLY : gauss_dist, set_random_seed 696 ! 697 IMPLICIT NONE 698 ! 699 INTEGER :: na, nb 700 REAL(DP) :: total_mass, kt, sigma, ek, ml(3), system_temp 701 ! 702 ! ... next command prevents different MD runs to start 703 ! ... with exactly the same "random" velocities 704 ! 705 CALL set_random_seed( ) 706 kt = temperature / ry_to_kelvin 707 ! 708 ! ... starting velocities have a Maxwell-Boltzmann distribution 709 ! 710 DO na = 1, nat 711 ! 712 sigma = SQRT( kt / mass(na) ) 713 ! 714 ! ... N.B. velocities must in a.u. units of alat 715 ! 716 vel(:,na) = gauss_dist( 0.D0, sigma, 3 ) / alat 717 ! 718 ENDDO 719 ! 720 ! ... the velocity of fixed ions must be zero 721 ! 722 vel = vel * DBLE( if_pos ) 723 ! 724 IF ( lconstrain ) THEN 725 ! 726 ! ... remove the component of the velocity along the 727 ! ... constraint gradient 728 ! 729 CALL remove_constr_vec( nat, tau, if_pos, ityp, alat, vel ) 730 ! 731 ENDIF 732 ! 733 IF ( invsym ) THEN 734 ! 735 ! ... if there is inversion symmetry, equivalent atoms have 736 ! ... opposite velocities 737 ! 738 DO na = 1, nat 739 ! 740 nb = irt( ( nsym / 2 + 1 ), na ) 741 ! 742 IF ( nb > na ) vel(:,nb) = - vel(:,na) 743 ! 744 ! ... the atom on the inversion center is kept fixed 745 ! 746 IF ( na == nb ) vel(:,na) = 0.D0 747 ! 748 ENDDO 749 ! 750 ELSE 751 ! 752 ! ... put total linear momentum equal zero if all atoms 753 ! ... are free to move 754 ! 755 ml(:) = 0.D0 756 ! 757 IF ( .NOT. ANY( if_pos(:,:) == 0 ) ) THEN 758 ! 759 total_mass = SUM( mass(1:nat) ) 760 DO na = 1, nat 761 ml(:) = ml(:) + mass(na)*vel(:,na) 762 ENDDO 763 ml(:) = ml(:) / total_mass 764 ! 765 ENDIF 766 ! 767 ENDIF 768 ! 769 ek = 0.D0 770 ! 771 DO na = 1, nat 772 ! 773 vel(:,na) = vel(:,na) - ml(:) 774 ! 775 ek = ek + 0.5D0 * mass(na) * & 776 ( ( vel(1,na) )**2 + ( vel(2,na) )**2 + ( vel(3,na) )**2 ) 777 ! 778 ENDDO 779 ! 780 ! ... after the velocity of the center of mass has been subtracted the 781 ! ... temperature is usually changed. Set again the temperature to the 782 ! ... right value. 783 ! 784 system_temp = 2.D0 / DBLE( ndof ) * ek * alat**2 * ry_to_kelvin 785 ! 786 CALL thermalize( 0, system_temp, temperature ) 787 ! 788 END SUBROUTINE start_therm 789 ! 790 END SUBROUTINE verlet 791 ! 792 ! 793 !------------------------------------------------------------------------ 794 SUBROUTINE terminate_verlet 795 !------------------------------------------------------------------------ 796 !! Terminate Verlet molecular dynamics calculation. 797 ! 798 USE io_global, ONLY : stdout 799 ! 800 WRITE( UNIT = stdout, & 801 FMT = '(/,5X,"The maximum number of steps has been reached.")' ) 802 WRITE( UNIT = stdout, & 803 FMT = '(/,5X,"End of molecular dynamics calculation")' ) 804 ! 805 CALL print_averages() 806 ! 807 END SUBROUTINE terminate_verlet 808 ! 809 ! 810 !------------------------------------------------------------------------ 811 SUBROUTINE proj_verlet( conv_ions ) 812 !------------------------------------------------------------------------ 813 !! This routine performs one step of structural relaxation using 814 !! the preconditioned-projected-Verlet algorithm. 815 ! 816 USE ions_base, ONLY : nat, ityp, tau, if_pos 817 USE cell_base, ONLY : alat 818 USE ener, ONLY : etot 819 USE force_mod, ONLY : force 820 USE relax, ONLY : epse, epsf 821 USE control_flags, ONLY : istep, lconstrain 822 ! 823 USE constraints_module, ONLY : remove_constr_force, check_constraint 824 ! TB 825 USE extfield, ONLY : relaxz 826 ! 827 IMPLICIT NONE 828 ! 829 LOGICAL, INTENT(OUT) :: conv_ions 830 ! 831 REAL(DP), ALLOCATABLE :: step(:,:) 832 REAL(DP) :: norm_step, etotold, delta(3) 833 INTEGER :: na 834 LOGICAL :: file_exists,leof 835 ! 836 REAL(DP), PARAMETER :: step_max = 0.6D0 ! bohr 837 ! 838 REAL(DP), EXTERNAL :: dnrm2 839 ! 840 ! 841 ALLOCATE( step( 3, nat ) ) 842 ! 843 tau_old(:,:) = tau(:,:) 844 tau_new(:,:) = 0.D0 845 vel(:,:) = 0.D0 846 acc(:,:) = 0.D0 847 conv_ions = .FALSE. 848 ! 849 CALL seqopn( 4, 'md', 'FORMATTED', file_exists ) 850 ! 851 IF ( file_exists ) THEN 852 ! 853 ! ... the file is read 854 ! 855 READ( UNIT = 4, FMT = * ) etotold, istep, tau_old(:,:) 856 ! 857 CLOSE( UNIT = 4, STATUS = 'KEEP' ) 858 ! 859 ELSE 860 ! 861 CLOSE( UNIT = 4, STATUS = 'DELETE' ) 862 ! 863 ! ... atoms are refold in the central box 864 ! 865 IF ( refold_pos ) CALL refold_tau() 866 ! 867 tau_old(:,:) = tau(:,:) 868 etotold = etot 869 istep = 0 870 WRITE( UNIT = stdout, & 871 FMT = '(/,5X,"Damped Dynamics Calculation")' ) 872 ! 873 ENDIF 874 ! 875 IF ( lconstrain ) THEN 876 ! 877 ! ... we first remove the component of the force along the 878 ! ... constraint gradient (this constitutes the initial guess 879 ! ... for the calculation of the lagrange multipliers) 880 ! 881 CALL remove_constr_force( nat, tau, if_pos, ityp, alat, force ) 882 ! 883#if ! defined (__REDUCE_OUTPUT) 884 ! 885 WRITE( stdout, '(/,5X,"Constrained forces (Ry/au):",/)') 886 ! 887 DO na = 1, nat 888 ! 889 WRITE( stdout, & 890 '(5X,"atom ",I3," type ",I2,3X,"force = ",3F14.8)' ) & 891 na, ityp(na), force(:,na) 892 ! 893 ENDDO 894 ! 895 WRITE( stdout, & 896 '(/5X,"Total force = ",F12.6)') dnrm2( 3*nat, force, 1 ) 897 ! 898#endif 899 ! 900 ENDIF 901 ! 902 ! ... check if convergence for structural minimization is achieved 903 ! 904 conv_ions = ( etotold - etot ) < epse 905 conv_ions = conv_ions .and. ( MAXVAL( ABS( force ) ) < epsf ) 906 ! 907 IF ( conv_ions ) THEN 908 ! 909 WRITE( UNIT = stdout, & 910 FMT = '(/,5X,"Damped Dynamics: convergence achieved in " & 911 & ,I3," steps")' ) istep 912 WRITE( UNIT = stdout, & 913 FMT = '(/,5X,"End of damped dynamics calculation")' ) 914 WRITE( UNIT = stdout, & 915 FMT = '(/,5X,"Final energy = ",F18.10," Ry"/)' ) etot 916 ! 917 CALL output_tau( .TRUE., .TRUE. ) 918 ! 919 RETURN 920 ! 921 ENDIF 922 ! 923 istep = istep + 1 924 WRITE( stdout, '(/,5X,"Entering Dynamics:",& 925 & T28,"iteration",T37," = ",I5)' ) istep 926 ! 927 ! ... Damped dynamics ( based on the projected-Verlet algorithm ) 928 ! 929 vel(:,:) = tau(:,:) - tau_old(:,:) 930 ! 931 CALL force_precond( istep, force, etotold ) 932 ! 933 acc(:,:) = force(:,:) / alat / amu_ry 934 ! 935 CALL project_velocity() 936 ! 937 step(:,:) = vel(:,:) + dt**2 * acc(:,:) 938 ! 939 norm_step = dnrm2( 3*nat, step, 1 ) 940 ! 941 step(:,:) = step(:,:) / norm_step 942 ! 943 tau_new(:,:) = tau(:,:) + step(:,:)*MIN( norm_step, step_max / alat ) 944 ! 945 ! TB 946 !IF ( .NOT. ANY( if_pos(:,:) == 0 ) ) THEN 947 IF ( .NOT. ANY( if_pos(:,:) == 0 ) .AND. (relaxz) ) THEN 948 WRITE( stdout, '("relaxz = .TRUE. => displacement of the center of mass is not subtracted")') 949 ENDIF 950 IF ( (.NOT. ANY( if_pos(:,:) == 0 )) .AND. (.NOT. relaxz) ) THEN 951 ! 952 ! ... if no atom has been fixed we compute the displacement of the 953 ! ... center of mass and we subtract it from the displaced positions 954 ! 955 delta(:) = 0.D0 956 ! 957 DO na = 1, nat 958 ! 959 delta(:) = delta(:) + ( tau_new(:,na) - tau(:,na) ) 960 ! 961 ENDDO 962 ! 963 delta(:) = delta(:) / DBLE( nat ) 964 ! 965 FORALL( na = 1:nat ) tau_new(:,na) = tau_new(:,na) - delta(:) 966 ! 967 ENDIF 968 ! 969 IF ( lconstrain ) THEN 970 ! 971 ! ... check if the new positions satisfy the constrain equation 972 ! 973 CALL check_constraint( nat, tau_new, tau, & 974 force, if_pos, ityp, alat, dt, amu_ry ) 975 ! 976 ENDIF 977 ! 978 ! ... save on file all the needed quantities 979 ! 980 CALL seqopn( 4, 'md', 'FORMATTED', file_exists ) 981 ! 982 leof = .TRUE. 983 WRITE( UNIT = 4, FMT = * ) etot, istep, tau(:,:), leof 984 ! 985 CLOSE( UNIT = 4, STATUS = 'KEEP' ) 986 ! 987 ! ... here the tau are shifted 988 ! 989 tau(:,:) = tau_new(:,:) 990 ! 991#if ! defined (__REDUCE_OUTPUT) 992 ! 993 CALL output_tau( .FALSE., .FALSE. ) 994 ! 995#endif 996 ! 997 DEALLOCATE( step ) 998 ! 999 END SUBROUTINE proj_verlet 1000 ! 1001 ! 1002 !------------------------------------------------------------------------ 1003 SUBROUTINE langevin_md() 1004 !------------------------------------------------------------------------ 1005 !! Langevin molecular dynamics. 1006 ! 1007 USE ions_base, ONLY : nat, ityp, tau, if_pos 1008 USE cell_base, ONLY : alat 1009 USE ener, ONLY : etot 1010 USE force_mod, ONLY : force 1011 USE control_flags, ONLY : istep, lconstrain 1012 USE random_numbers, ONLY : gauss_dist 1013 ! 1014 USE constraints_module, ONLY : nconstr 1015 USE constraints_module, ONLY : remove_constr_force, check_constraint 1016 ! 1017 IMPLICIT NONE 1018 ! 1019 REAL(DP) :: sigma, kt 1020 REAL(DP) :: delta(3) 1021 INTEGER :: na 1022 LOGICAL :: file_exists 1023 ! 1024 REAL(DP), EXTERNAL :: dnrm2 1025 ! 1026 CALL seqopn( 4, 'md', 'FORMATTED', file_exists ) 1027 ! 1028 IF ( file_exists ) THEN 1029 ! 1030 ! ... the file is read : simulation is continuing 1031 ! 1032 READ( UNIT = 4, FMT = * ) istep 1033 ! 1034 CLOSE( UNIT = 4, STATUS = 'KEEP' ) 1035 ! 1036 ELSE 1037 ! 1038 CLOSE( UNIT = 4, STATUS = 'DELETE' ) 1039 ! 1040 ! ... the file is absent : simulation is starting from scratch 1041 ! 1042 istep = 0 1043 ! 1044 WRITE( UNIT = stdout, & 1045 FMT = '(/,5X,"Over-damped Langevin Dynamics Calculation")' ) 1046 ! 1047 ! ... atoms are refold in the central box if required 1048 ! 1049 IF ( refold_pos ) CALL refold_tau() 1050 ! 1051 WRITE( UNIT = stdout, & 1052 FMT = '(5X,"Integration step",T27," = ",F8.2," a.u.,")' ) dt 1053 ! 1054 ENDIF 1055 ! 1056 istep = istep + 1 1057 WRITE( UNIT = stdout, & 1058 FMT = '(/,5X,"Entering Dynamics:",T28, & 1059 & "iteration",T37," = ",I5,/)' ) istep 1060 ! 1061 IF ( lconstrain ) THEN 1062 ! 1063 ! ... we first remove the component of the force along the 1064 ! ... constraint gradient ( this constitutes the initial 1065 ! ... guess for the calculation of the lagrange multipliers ) 1066 ! 1067 CALL remove_constr_force( nat, tau, if_pos, ityp, alat, force ) 1068 ! 1069 ENDIF 1070 ! 1071 ! ... compute the stochastic term 1072 ! 1073 kt = temperature / ry_to_kelvin 1074 ! 1075 sigma = SQRT( 2.D0*dt*kt ) 1076 ! 1077 delta(:) = 0.D0 1078 ! 1079 DO na = 1, nat 1080 ! 1081 chi(:,na) = gauss_dist( 0.D0, sigma, 3 )*DBLE( if_pos(:,na) ) 1082 ! 1083 delta(:) = delta(:) + chi(:,na) 1084 ! 1085 ENDDO 1086 ! 1087 FORALL( na = 1:nat ) chi(:,na) = chi(:,na) - delta(:) / DBLE( nat ) 1088 ! 1089 PRINT *, "|F| = ", dt*dnrm2( 3*nat, force, 1 ) 1090 PRINT *, "|CHI| = ", dnrm2( 3*nat, chi, 1 ) 1091 ! 1092 ! ... over-damped Langevin dynamics 1093 ! 1094 tau_new(:,:) = tau(:,:) + ( dt*force(:,:) + chi(:,:) ) / alat 1095 ! 1096 IF ( .NOT. ANY( if_pos(:,:) == 0 ) ) THEN 1097 ! 1098 ! ... here we compute the displacement of the center of mass and we 1099 ! ... subtract it from the displaced positions 1100 ! 1101 delta(:) = 0.D0 1102 ! 1103 DO na = 1, nat 1104 ! 1105 delta(:) = delta(:) + ( tau_new(:,na) - tau(:,na) ) 1106 ! 1107 ENDDO 1108 ! 1109 FORALL( na = 1:nat ) tau_new(:,na) = tau_new(:,na) - delta(:) 1110 ! 1111 ENDIF 1112 ! 1113 IF ( lconstrain ) THEN 1114 ! 1115 ! ... check if the new positions satisfy the constrain equation 1116 ! 1117 CALL check_constraint( nat, tau_new, tau, & 1118 force, if_pos, ityp, alat, dt, amu_ry ) 1119 ! 1120#if ! defined (__REDUCE_OUTPUT) 1121 ! 1122 WRITE( stdout, '(/,5X,"Constrained forces (Ry/au):",/)') 1123 ! 1124 DO na = 1, nat 1125 ! 1126 WRITE( stdout, & 1127 '(5X,"atom ",I3," type ",I2,3X,"force = ",3F14.8)' ) & 1128 na, ityp(na), force(:,na) 1129 ! 1130 ENDDO 1131 ! 1132 WRITE( stdout, '(/5X,"Total force = ",F12.6)') dnrm2( 3*nat, force, 1 ) 1133 ! 1134#endif 1135 ! 1136 ENDIF 1137 ! 1138 ! ... save all the needed quantities on file 1139 ! 1140 CALL seqopn( 4, 'md', 'FORMATTED', file_exists ) 1141 ! 1142 WRITE( UNIT = 4, FMT = * ) istep 1143 ! 1144 CLOSE( UNIT = 4, STATUS = 'KEEP' ) 1145 ! 1146 ! ... here the tau are shifted 1147 ! 1148 tau(:,:) = tau_new(:,:) 1149 ! 1150#if ! defined (__REDUCE_OUTPUT) 1151 ! 1152 CALL output_tau( .FALSE., .FALSE. ) 1153 ! 1154#endif 1155 ! 1156 END SUBROUTINE langevin_md 1157 ! 1158 ! 1159 !----------------------------------------------------------------------- 1160 SUBROUTINE refold_tau() 1161 !----------------------------------------------------------------------- 1162 !! Refold atomic positions. 1163 ! 1164 USE ions_base, ONLY : nat, tau 1165 USE cell_base, ONLY : alat 1166 USE constraints_module, ONLY : pbc 1167 ! 1168 IMPLICIT NONE 1169 ! 1170 INTEGER :: ia 1171 ! 1172 ! 1173 DO ia = 1, nat 1174 ! 1175 tau(:,ia) = pbc( tau(:,ia) * alat ) / alat 1176 ! 1177 ENDDO 1178 ! 1179 END SUBROUTINE refold_tau 1180 ! 1181 ! 1182 !----------------------------------------------------------------------- 1183 SUBROUTINE compute_averages( istep ) 1184 !----------------------------------------------------------------------- 1185 !! Molecular dynamics - compute averages. 1186 ! 1187 USE ions_base, ONLY : nat, tau, fixatom 1188 USE cell_base, ONLY : alat, at 1189 USE constraints_module, ONLY : pbc 1190 USE io_files, ONLY : delete_if_present 1191 ! 1192 IMPLICIT NONE 1193 ! 1194 INTEGER, INTENT(in) :: istep 1195 !! md step 1196 ! 1197 ! ... local variables 1198 ! 1199 INTEGER :: i, j, idx 1200 REAL(DP) :: dx, dy, dz 1201 REAL(DP) :: dtau(3) 1202 REAL(DP) :: inv_dmax 1203 REAL(DP), ALLOCATABLE :: msd(:) 1204 REAL(DP), PARAMETER :: max_dist(3) = (/ 0.5D0, 0.5D0, 0.5D0 /) 1205 ! 1206 ! ... MSD and diffusion coefficient 1207 ! 1208 ALLOCATE( msd( nat ) ) 1209 ! 1210 IF ( istep == 1 ) THEN 1211 ! 1212 radial_distr(:,:) = 0.D0 1213 ! 1214 CALL delete_if_present( TRIM( tmp_dir ) // & 1215 & TRIM( prefix ) // ".msd.dat" ) 1216 ! 1217 ENDIF 1218 ! 1219 DO i = 1, nat 1220 ! 1221 dx = ( tau(1,i) - tau_ref(1,i) ) * alat 1222 dy = ( tau(2,i) - tau_ref(2,i) ) * alat 1223 dz = ( tau(3,i) - tau_ref(3,i) ) * alat 1224 ! 1225 msd(i) = dx*dx + dy*dy + dz*dz 1226 ! 1227 ENDDO 1228 ! 1229 diff_coeff(:) = msd(:) / ( 6.D0*DBLE( istep )*dt ) 1230 ! 1231 ! ... conversion from Rydberg atomic units to cm^2/sec 1232 ! 1233 diff_coeff(:) = diff_coeff(:) * bohr_radius_cm**2 / ( 2.D-12*au_ps ) 1234 ! 1235 OPEN( UNIT = 4, POSITION = 'APPEND', & 1236 FILE = TRIM( tmp_dir ) // TRIM( prefix ) // ".msd.dat" ) 1237 ! 1238 WRITE( 4, '(2(2X,F16.8))' ) & 1239 ( istep*dt*2.D0*au_ps ), SUM( msd(:) ) / DBLE( nat-fixatom ) 1240 ! 1241 CLOSE( UNIT = 4, STATUS = 'KEEP' ) 1242 ! 1243 DEALLOCATE( msd ) 1244 ! 1245 ! ... radial distribution function g(r) 1246 ! 1247 inv_dmax = 1.D0 / ( norm( MATMUL( at(:,:), max_dist(:) ) ) * alat ) 1248 ! 1249 DO i = 1, nat 1250 ! 1251 DO j = 1, nat 1252 ! 1253 IF ( i == j ) CYCLE 1254 ! 1255 dtau(:) = pbc( ( tau(:,i) - tau(:,j) ) * alat ) 1256 ! 1257 idx = ANINT( norm( dtau(:) ) * inv_dmax * DBLE( hist_len ) ) 1258 ! 1259 IF( idx > 0 .and. idx <= SIZE( radial_distr, 1 ) ) & 1260 radial_distr(idx,i) = radial_distr(idx,i) + 1.D0 1261 ! 1262 ENDDO 1263 ! 1264 ENDDO 1265 ! 1266 END SUBROUTINE compute_averages 1267 ! 1268 ! 1269 !----------------------------------------------------------------------- 1270 SUBROUTINE print_averages() 1271 !----------------------------------------------------------------------- 1272 !! Molecular dynamics - print averages. 1273 ! 1274 USE control_flags, ONLY : nstep 1275 USE cell_base, ONLY : omega, at, alat 1276 USE ions_base, ONLY : nat, fixatom 1277 ! 1278 IMPLICIT NONE 1279 ! 1280 INTEGER :: i, idx 1281 REAL(DP) :: dist, dmax 1282 REAL(DP), PARAMETER :: max_dist(3) = (/ 0.5D0, 0.5D0, 0.5D0 /) 1283 ! 1284 ! ... diffusion coefficient 1285 ! 1286 WRITE( UNIT = stdout, & 1287 FMT = '(/,5X,"diffusion coefficients :")' ) 1288 ! 1289 DO i = 1, nat 1290 ! 1291 WRITE( UNIT = stdout, & 1292 FMT = '(5X,"atom ",I5," D = ",F16.8," cm^2/s")' ) & 1293 i, diff_coeff(i) 1294 ! 1295 ENDDO 1296 ! 1297 WRITE( UNIT = stdout, FMT = '(/,5X,"< D > = ",F16.8," cm^2/s")' ) & 1298 sum( diff_coeff(:) ) / DBLE( nat-fixatom ) 1299 ! 1300 ! ... radial distribution function g(r) 1301 ! 1302 dmax = norm( matmul( at(:,:), max_dist(:) ) ) * alat 1303 ! 1304 radial_distr(:,:) = radial_distr(:,:) * omega / DBLE( nat ) / fpi 1305 ! 1306 radial_distr(:,:) = radial_distr(:,:) / ( dmax / DBLE( hist_len ) ) 1307 ! 1308 radial_distr(:,:) = radial_distr(:,:) / DBLE( nstep ) 1309 ! 1310 OPEN( UNIT = 4, FILE = TRIM( tmp_dir ) // TRIM( prefix ) // ".rdf.dat" ) 1311 ! 1312 DO idx = 1, hist_len 1313 ! 1314 dist = DBLE( idx ) / DBLE( hist_len ) * dmax 1315 ! 1316 IF ( dist > dmax / SQRT( 3.0d0 ) ) CYCLE 1317 ! 1318 radial_distr(idx,:) = radial_distr(idx,:) / dist**2 1319 ! 1320 WRITE( 4, '(2(2X,F16.8))' ) & 1321 dist, SUM( radial_distr(idx,:) ) / DBLE( nat ) 1322 ! 1323 ENDDO 1324 ! 1325 CLOSE( UNIT = 4 ) 1326 ! 1327 END SUBROUTINE print_averages 1328 ! 1329 ! 1330 !----------------------------------------------------------------------- 1331 SUBROUTINE force_precond( istep, force, etotold ) 1332 !----------------------------------------------------------------------- 1333 !! This routine computes an estimate of \(H^{-1}\) by using the BFGS 1334 !! algorithm and the preconditioned gradient \(\text{pg} = H^{-1} g\). 1335 !! It works in atomic units. 1336 ! 1337 USE ener, ONLY : etot 1338 USE cell_base, ONLY : alat 1339 USE ions_base, ONLY : nat, tau 1340 USE io_files, ONLY : iunbfgs, tmp_dir 1341 ! 1342 IMPLICIT NONE 1343 ! 1344 INTEGER, INTENT(IN) :: istep 1345 REAL(DP), INTENT(INOUT) :: force(:,:) 1346 REAL(DP), INTENT(IN) :: etotold 1347 ! 1348 REAL(DP), ALLOCATABLE :: pos(:), pos_p(:) 1349 REAL(DP), ALLOCATABLE :: grad(:), grad_p(:), precond_grad(:) 1350 REAL(DP), ALLOCATABLE :: inv_hess(:,:) 1351 REAL(DP), ALLOCATABLE :: y(:), s(:) 1352 REAL(DP), ALLOCATABLE :: Hy(:), yH(:) 1353 REAL(DP) :: sdoty, pg_norm 1354 INTEGER :: dim 1355 CHARACTER(LEN=256) :: bfgs_file 1356 LOGICAL :: file_exists 1357 ! 1358 INTEGER, PARAMETER :: nrefresh = 25 1359 REAL(DP), PARAMETER :: max_pg_norm = 0.8D0 1360 ! 1361 ! 1362 dim = 3 * nat 1363 ! 1364 ALLOCATE( pos( dim ), pos_p( dim ) ) 1365 ALLOCATE( grad( dim ), grad_p( dim ), precond_grad( dim ) ) 1366 ALLOCATE( y( dim ), s( dim ) ) 1367 ALLOCATE( inv_hess( dim, dim ) ) 1368 ALLOCATE( Hy( dim ), yH( dim ) ) 1369 ! 1370 pos(:) = RESHAPE( tau, (/ dim /) ) * alat 1371 grad(:) = - RESHAPE( force, (/ dim /) ) 1372 ! 1373 bfgs_file = TRIM( tmp_dir ) // TRIM( prefix ) // '.bfgs' 1374 ! 1375 INQUIRE( FILE = TRIM( bfgs_file ) , EXIST = file_exists ) 1376 ! 1377 IF ( file_exists ) THEN 1378 ! 1379 OPEN( UNIT = iunbfgs, & 1380 FILE = TRIM( bfgs_file ), STATUS = 'OLD', ACTION = 'READ' ) 1381 ! 1382 READ( iunbfgs, * ) pos_p 1383 READ( iunbfgs, * ) grad_p 1384 READ( iunbfgs, * ) inv_hess 1385 ! 1386 CLOSE( UNIT = iunbfgs ) 1387 ! 1388 ! ... the approximate inverse hessian is reset to one every nrefresh 1389 ! ... iterations: this is one to clean-up the memory of the starting 1390 ! ... configuration 1391 ! 1392 IF ( MOD( istep, nrefresh ) == 0 ) inv_hess(:,:) = identity( dim ) 1393 ! 1394 IF ( etot < etotold ) THEN 1395 ! 1396 ! ... BFGS update 1397 ! 1398 s(:) = pos(:) - pos_p(:) 1399 y(:) = grad(:) - grad_p(:) 1400 ! 1401 sdoty = ( s(:) .DOT. y(:) ) 1402 ! 1403 IF ( sdoty > eps8 ) THEN 1404 ! 1405 Hy(:) = ( inv_hess(:,:) .TIMES. y(:) ) 1406 yH(:) = ( y(:) .TIMES. inv_hess(:,:) ) 1407 ! 1408 inv_hess = inv_hess + 1.D0 / sdoty * & 1409 ( ( 1.D0 + ( y .DOT. Hy ) / sdoty ) * matrix( s, s ) - & 1410 ( matrix( s, yH ) + matrix( Hy, s ) ) ) 1411 ! 1412 ENDIF 1413 ! 1414 ENDIF 1415 ! 1416 ELSE 1417 ! 1418 inv_hess(:,:) = identity( dim ) 1419 ! 1420 ENDIF 1421 ! 1422 precond_grad(:) = ( inv_hess(:,:) .TIMES. grad(:) ) 1423 ! 1424 IF ( ( precond_grad(:) .DOT. grad(:) ) < 0.D0 ) THEN 1425 ! 1426 WRITE( UNIT = stdout, & 1427 FMT = '(/,5X,"uphill step: resetting bfgs history",/)' ) 1428 ! 1429 precond_grad(:) = grad(:) 1430 ! 1431 inv_hess(:,:) = identity( dim ) 1432 ! 1433 ENDIF 1434 ! 1435 OPEN( UNIT = iunbfgs, & 1436 FILE = TRIM( bfgs_file ), STATUS = 'UNKNOWN', ACTION = 'WRITE' ) 1437 ! 1438 WRITE( iunbfgs, * ) pos(:) 1439 WRITE( iunbfgs, * ) grad(:) 1440 WRITE( iunbfgs, * ) inv_hess(:,:) 1441 ! 1442 CLOSE( UNIT = iunbfgs ) 1443 ! 1444 ! ... the length of the step is always shorter than pg_norm 1445 ! 1446 pg_norm = norm( precond_grad(:) ) 1447 ! 1448 precond_grad(:) = precond_grad(:) / pg_norm 1449 precond_grad(:) = precond_grad(:) * MIN( pg_norm, max_pg_norm ) 1450 ! 1451 force(:,:) = - RESHAPE( precond_grad(:), (/ 3, nat /) ) 1452 ! 1453 DEALLOCATE( pos, pos_p ) 1454 DEALLOCATE( grad, grad_p, precond_grad ) 1455 DEALLOCATE( inv_hess ) 1456 DEALLOCATE( y, s ) 1457 DEALLOCATE( Hy, yH ) 1458 ! 1459 END SUBROUTINE force_precond 1460 ! 1461 !----------------------------------------------------------------------- 1462 SUBROUTINE project_velocity() 1463 !----------------------------------------------------------------------- 1464 !! Quick-min algorithm. 1465 ! 1466 USE control_flags, ONLY : istep 1467 USE ions_base, ONLY : nat 1468 ! 1469 IMPLICIT NONE 1470 ! 1471 REAL(DP) :: norm_acc, projection 1472 REAL(DP), ALLOCATABLE :: acc_versor(:,:) 1473 ! 1474 REAL(DP), EXTERNAL :: dnrm2, ddot 1475 ! 1476 ! 1477 IF ( istep == 1 ) RETURN 1478 ! 1479 ALLOCATE( acc_versor( 3, nat ) ) 1480 ! 1481 norm_acc = dnrm2( 3*nat, acc(:,:), 1 ) 1482 ! 1483 acc_versor(:,:) = acc(:,:) / norm_acc 1484 ! 1485 projection = ddot( 3*nat, vel(:,:), 1, acc_versor(:,:), 1 ) 1486 ! 1487 WRITE( UNIT = stdout, FMT = '(/,5X,"<vel(dt)|acc(dt)> = ",F12.8)' ) & 1488 projection / dnrm2( 3*nat, vel, 1 ) 1489 ! 1490 vel(:,:) = acc_versor(:,:) * MAX( 0.D0, projection ) 1491 ! 1492 DEALLOCATE( acc_versor ) 1493 ! 1494 END SUBROUTINE project_velocity 1495 ! 1496 !----------------------------------------------------------------------- 1497 SUBROUTINE thermalize( nraise, system_temp, required_temp ) 1498 !----------------------------------------------------------------------- 1499 !! * Berendsen rescaling (Eq. 7.59 of Allen & Tildesley); 1500 !! * rescale the velocities by a factor 3 / 2KT / Ek. 1501 ! 1502 IMPLICIT NONE 1503 ! 1504 REAL(DP), INTENT(IN) :: system_temp, required_temp 1505 INTEGER, INTENT(IN) :: nraise 1506 ! 1507 REAL(DP) :: aux 1508 ! 1509 IF ( nraise > 0 ) THEN 1510 ! 1511 ! ... Berendsen rescaling (Eq. 7.59 of Allen & Tildesley) 1512 ! ... the "rise time" is tau=nraise*dt so dt/tau=1/nraise 1513 ! ... Equivalent to traditional rescaling if nraise=1 1514 ! 1515 IF ( system_temp > 0.D0 .AND. required_temp > 0.D0 ) THEN 1516 aux = SQRT( 1.d0 + (required_temp / system_temp - 1.d0) * & 1517 (1.D0/DBLE(nraise) ) ) 1518 ELSE 1519 aux = 0.d0 1520 ENDIF 1521 ! 1522 ELSE 1523 ! 1524 ! ... rescale the velocities by a factor 3 / 2KT / Ek 1525 ! 1526 IF ( system_temp > 0.D0 .AND. required_temp > 0.D0 ) THEN 1527 aux = SQRT( required_temp / system_temp ) 1528 ELSE 1529 aux = 0.d0 1530 ENDIF 1531 ! 1532 ENDIF 1533 ! 1534 vel(:,:) = vel(:,:) * aux 1535 ! 1536 END SUBROUTINE thermalize 1537 ! 1538 ! 1539 !----------------------------------------------------------------------- 1540 SUBROUTINE smart_MC() 1541 !----------------------------------------------------------------------- 1542 !! Routine to apply smart_MC. 1543 !! Implemented by Xiaochuan Ge, Jul., 2013 1544 ! 1545 !! At this moment works only with Langevin dynamics. 1546 !! For the formula see R.J.Rossky, JCP, 69, 4628(1978). 1547 ! 1548 USE ions_base, ONLY : nat, ityp, tau, if_pos,atm 1549 USE cell_base, ONLY : alat 1550 USE ener, ONLY : etot 1551 USE force_mod, ONLY : force 1552 USE control_flags, ONLY : istep, lconstrain 1553 USE constraints_module, ONLY : remove_constr_force, check_constraint 1554 USE random_numbers, ONLY : randy 1555 USE io_files, ONLY : prefix 1556 USE io_global, ONLY : ionode 1557 USE constants, ONLY : bohr_radius_angs 1558 ! 1559 IMPLICIT NONE 1560 ! 1561 LOGICAL :: accept 1562 REAL(DP) :: kt,sigma2, & 1563 T_ij,T_ji,boltzman_ji, & ! boltzman_ji=exp[-(etot_new-etot_old)/kt] 1564 temp,p_smc ! *_smart means *_old, the quantity of the 1565 ! previous step 1566 ! 1567 INTEGER :: ia, ip 1568 ! 1569 IF ( lconstrain ) THEN 1570 ! ... we first remove the component of the force along the 1571 ! ... constraint gradient ( this constitutes the initial 1572 ! ... guess for the calculation of the lagrange multipliers ) 1573 CALL remove_constr_force( nat, tau, if_pos, ityp, alat, force ) 1574 ENDIF 1575 ! 1576 IF (first_iter) THEN ! For the first iteration 1577 ALLOCATE( tau_smart(3,nat) ) 1578 ALLOCATE( force_smart(3,nat) ) 1579 tau_smart = tau 1580 etot_smart = etot 1581 force_smart = force 1582 first_iter = .FALSE. 1583 RETURN 1584 ENDIF 1585 ! 1586 kt = temperature / ry_to_kelvin 1587 sigma2 = 2.D0*dt*kt 1588 ! 1589 T_ij=0.0d0 1590 T_ji=0.0d0 1591 DO ia = 1, nat 1592 DO ip = 1, 3 1593 T_ij = T_ij + ( (tau(ip,ia)-tau_smart(ip,ia))*alat-dt*force_smart(ip,ia) )**2 1594 T_ji = T_ji + ( (tau_smart(ip,ia)-tau(ip,ia))*alat-dt*force(ip,ia) )**2 1595 ENDDO 1596 ENDDO 1597 T_ij = EXP(-T_ij/(2*sigma2)) 1598 T_ji = EXP(-T_ji/(2*sigma2)) 1599 ! 1600 boltzman_ji = EXP(-(etot-etot_smart)/kt) 1601 ! 1602 p_smc = T_ji*boltzman_ji/T_ij 1603 ! 1604 WRITE(stdout, '(5x,"The old energy is:",3x,F17.8," Ry")') etot_smart 1605 WRITE(stdout, '(5x,"The new energy is:",3x,F17.8," Ry")') etot 1606 WRITE(stdout, '(5x,"The possibility to accept this step is:",3x,F10.7/)') p_smc 1607 WRITE(stdout, '(5x,"Nervously waiting for the fate ..."/)') 1608 ! 1609 ! Decide if accept the new config 1610 temp = randy() 1611 WRITE(stdout, '(5x,"The fate says:",5x,F10.7)') temp 1612 IF(temp <= p_smc) THEN 1613 WRITE(stdout, '(5x,"The new config is accepted")') 1614 num_accept=num_accept+1 1615 tau_smart=tau 1616 etot_smart=etot 1617 force_smart=force 1618 ELSE 1619 WRITE(stdout, '(5x,"The new config is not accepted")') 1620 tau=tau_smart 1621 etot=etot_smart 1622 force=force_smart 1623 ENDIF 1624 ! 1625 WRITE (stdout, '(5x,"The current acceptance is :",3x,F10.6)') DBLE(num_accept)/istep 1626 ! 1627 ! Print the trajectory 1628#if defined(__MPI) 1629 IF(ionode) THEN 1630#endif 1631 OPEN(117,file="trajectory-"//TRIM(prefix)//".xyz", STATUS="unknown", POSITION='APPEND') 1632 WRITE(117,'(I5)') nat 1633 WRITE(117,'("# Step: ",I5,5x,"Total energy: ",F17.8,5x,"Ry")') istep-1, etot 1634 DO ia = 1, nat 1635 WRITE( 117, '(A3,3X,3F14.9)') atm(ityp(ia)),tau(:,ia)*alat*bohr_radius_angs 1636 ENDDO 1637 CLOSE(117) 1638#if defined(__MPI) 1639 ENDIF 1640#endif 1641 ! 1642 RETURN 1643 ! 1644 END SUBROUTINE smart_MC 1645 ! 1646 !----------------------------------------------------------------------- 1647 SUBROUTINE thermalize_resamp_vscaling( nraise, system_temp, required_temp ) 1648 !----------------------------------------------------------------------- 1649 !! Sample velocities using stochastic velocity rescaling, based on: 1650 !! Bussi, Donadio, Parrinello, J. Chem. Phys. 126, 014101 (2007), 1651 !! doi: 10.1063/1.2408420 1652 ! 1653 !! Implemented (2019) by Leonid Kahle and Ngoc Linh Nguyen, 1654 !! Theory and Simulations of Materials Laboratory, EPFL. 1655 ! 1656 USE ions_base, ONLY : nat, if_pos 1657 USE constraints_module, ONLY : nconstr 1658 USE cell_base, ONLY : alat 1659 USE random_numbers, ONLY : gauss_dist, sum_of_gaussians2 1660 ! 1661 IMPLICIT NONE 1662 ! 1663 REAL(DP), INTENT(in) :: system_temp, required_temp 1664 INTEGER, INTENT(in) :: nraise 1665 ! 1666 INTEGER :: i, ndof 1667 REAL(DP) :: factor, rr 1668 REAL(DP) :: aux, aux2 1669 real(DP), external :: gasdev, sumnoises 1670 INTEGER :: na 1671 ! 1672 ! ... the number of degrees of freedom 1673 ! 1674 IF ( ANY( if_pos(:,:) == 0 ) ) THEN 1675 ! 1676 ndof = 3*nat - count( if_pos(:,:) == 0 ) - nconstr 1677 ! 1678 ELSE 1679 ! 1680 ndof = 3*nat - 3 - nconstr 1681 ! 1682 ENDIF 1683 ! 1684 IF ( nraise > 0 ) THEN 1685 ! 1686 ! ... the "rise time" is tau=nraise*dt so dt/tau=1/nraise 1687 ! ... Equivalent to traditional rescaling if nraise=1 1688 ! 1689 factor = exp(-1.0/nraise) 1690 ELSE 1691 ! 1692 factor = 0.0 1693 ! 1694 ENDIF 1695 ! 1696 IF ( system_temp > 0.D0 .and. required_temp > 0.D0 ) THEN 1697 ! 1698 ! Applying Eq. (A7) from J. Chem. Phys. 126, 014101 (2007) 1699 ! 1700 rr = gauss_dist(0.0D0, 1.0D0) 1701 aux2 = factor + (1.0D0-factor)*( sum_of_gaussians2(ndof-1) +rr**2) & 1702 * required_temp/(ndof*system_temp) & 1703 + 2*rr*sqrt((factor*(1.0D0-factor)*required_temp)/(ndof*system_temp)) 1704 ! 1705 aux = sqrt(aux2) 1706 1707 ELSE 1708 ! 1709 aux = 0.d0 1710 ! 1711 ENDIF 1712 ! 1713 ! Global rescaling applied to velocities 1714 vel(:,:) = vel(:,:) * aux 1715 ! 1716 END SUBROUTINE thermalize_resamp_vscaling 1717 ! 1718END MODULE dynamics_module 1719