1! 2! Copyright (C) 2003-2007 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!--------------------------------------------------------------------------- 10MODULE path_base 11 !--------------------------------------------------------------------------- 12 ! 13 ! ... This module contains most of the subroutines and functions needed by 14 ! ... the implementation of "NEB" and "SMD" methods into Quantum ESPRESSO 15 ! 16 ! ... Other relevant files are: 17 ! 18 ! ... path_variables.f90 19 ! ... path_io_routines.f90 20 ! ... path_opt_routines.f90 21 ! ... path_reparametrisation.f90 22 ! ... path_formats.f90 23 ! ... compute_scf.f90 24 ! 25 ! ... The code is based on the NEB algorithm described in : 26 ! 27 ! ... 1) G. Henkelman, B.P. Uberuaga, and H. Jonsson; 28 ! ... J.Chem.Phys., 113, 9901, (2000) 29 ! ... 2) G. Henkelman, and H. Jonsson; J.Chem.Phys., 113, 9978, (2000) 30 ! 31 ! ... More details about the implementation can be found at 32 ! 33 ! ... http://www.sissa.it/cm/thesis/2005/sbraccia.pdf 34 ! 35 ! ... Code written and maintained by Carlo Sbraccia ( 2003-2007 ) 36 ! 37 USE kinds, ONLY : DP 38 USE constants, ONLY : eps32, pi, autoev, bohr_radius_angs, eV_to_kelvin 39 USE path_io_units_module, ONLY : iunpath 40 USE io_global, ONLY : meta_ionode, meta_ionode_id 41 USE mp, ONLY : mp_bcast 42 USE mp_world, ONLY : world_comm 43 ! 44 USE basic_algebra_routines 45 ! 46 PRIVATE 47 ! 48 PUBLIC :: initialize_path 49 PUBLIC :: search_mep 50 ! 51 CONTAINS 52 ! 53 ! ... module procedures 54 ! 55 !----------------------------------------------------------------------- 56 SUBROUTINE initialize_path() 57 !----------------------------------------------------------------------- 58 ! 59 USE control_flags, ONLY : conv_elec 60 USE ions_base, ONLY : amass, ityp 61 USE io_files, ONLY : prefix, tmp_dir 62 USE mp_images, ONLY : nimage 63 USE path_input_parameters_module, ONLY : pos_ => pos, & 64 climbing_ => climbing, & 65 input_images, nstep_path_ => nstep_path 66 USE path_input_parameters_module, ONLY : restart_mode 67 USE path_input_parameters_module, ONLY : nat 68 USE path_variables, ONLY : fix_atom_pos 69 USE path_variables, ONLY : climbing, pos, istep_path, nstep_path, & 70 dim1, num_of_images, pes, grad_pes, mass, & 71 use_masses, tangent, error, path_length, & 72 deg_of_freedom, frozen, use_freezing, k, & 73 k_min, tune_load_balance, grad, posold, & 74 elastic_grad, pending_image, first_last_opt 75 USE path_variables, ONLY : path_allocation 76 USE path_io_routines, ONLY : read_restart 77 USE path_io_units_module, ONLY : path_file, dat_file, crd_file, & 78 int_file, xyz_file, axsf_file, broy_file 79 USE fcp_variables, ONLY : lfcpopt 80 USE fcp_opt_routines, ONLY : fcp_opt_allocation 81 ! 82 IMPLICIT NONE 83 ! 84 INTEGER :: i, fii, lii 85 LOGICAL :: file_exists 86 ! 87 ! ... output files are set 88 ! 89 path_file = TRIM( prefix ) // ".path" 90 dat_file = TRIM( prefix ) // ".dat" 91 int_file = TRIM( prefix ) // ".int" 92 crd_file = TRIM( prefix ) // ".crd" 93 xyz_file = TRIM( prefix ) // ".xyz" 94 axsf_file = TRIM( prefix ) // ".axsf" 95 ! 96 broy_file = TRIM( tmp_dir ) // TRIM( prefix ) // ".broyden" 97 ! 98 ! ... istep_path is initialised to zero 99 ! 100 istep_path = 0 101 pending_image = 0 102 conv_elec = .TRUE. 103 ! 104 ! ... the dimension of all "path" arrays (dim1) is set here 105 ! ... ( it corresponds to the dimension of the configurational space ) 106 ! 107 ! 108 dim1 = 3*nat 109 ! 110 ! 111 IF ( nimage > 1 ) THEN 112 ! 113 ! ... the automatic tuning of the load balance in 114 ! ... image-parallelisation is switched off by default 115 ! 116 tune_load_balance = .FALSE. 117 ! 118 ! ... in the case of image-parallelisation the number of images 119 ! ... to be optimised must be larger than nimage 120 ! 121 IF ( first_last_opt ) THEN 122 ! 123 fii = 1 124 lii = num_of_images 125 ! 126 ELSE 127 ! 128 fii = 2 129 lii = num_of_images - 1 130 ! 131 END IF 132 ! 133 IF ( nimage > ( lii - fii + 1 ) ) & 134 CALL errore( 'initialize_path', 'nimage is ' // & 135 & 'larger than the available number of images', 1 ) 136 ! 137 END IF 138 ! 139 ! ... dynamical allocation of arrays 140 ! 141 CALL path_allocation() 142 if ( lfcpopt ) CALL fcp_opt_allocation() 143 ! 144 IF ( use_masses ) THEN 145 ! 146 ! ... mass weighted coordinates are used 147 ! 148 DO i = 1, nat 149 ! 150 mass(3*i-2) = amass(ityp(i)) 151 mass(3*i-1) = amass(ityp(i)) 152 mass(3*i-0) = amass(ityp(i)) 153 ! 154 END DO 155 ! 156 ELSE 157 ! 158 mass = 1.0_DP 159 ! 160 END IF 161 ! 162 ! ... initialization of the allocatable arrays 163 ! 164 pos(:,1:input_images) = pos_(1:dim1,1:input_images) 165 ! 166 pes = 0.0_DP 167 grad_pes = 0.0_DP 168 elastic_grad = 0.0_DP 169 tangent = 0.0_DP 170 grad = 0.0_DP 171 error = 0.0_DP 172 frozen = .FALSE. 173 ! 174 k = k_min 175 ! 176 IF ( ALLOCATED( climbing_ ) ) THEN 177 ! 178 climbing = climbing_ 179 ! 180 ELSE 181 ! 182 climbing = .FALSE. 183 ! 184 END IF 185 ! 186 ! ... initial path is read from file ( restart_mode == "restart" ) or 187 ! ... generated from the input images ( restart_mode = "from_scratch" ) 188 ! ... It is always read from file in the case of "free-energy" 189 ! ... calculations 190 ! 191 IF ( restart_mode == "restart" ) THEN 192 ! 193 IF ( meta_ionode ) THEN 194 ! 195 INQUIRE( FILE = path_file, EXIST = file_exists ) 196 ! 197 IF ( .NOT. file_exists ) THEN 198 ! 199 WRITE( iunpath, & 200 & '(/,5X,"restart file ''",A,"'' not found: ", & 201 & /,5X,"starting from scratch")' ) TRIM( path_file ) 202 ! 203 restart_mode = "from_scratch" 204 ! 205 END IF 206 ! 207 END IF 208 ! 209 CALL mp_bcast( restart_mode, meta_ionode_id, world_comm ) 210 ! 211 END IF 212 ! 213 IF ( restart_mode == "restart" ) THEN 214 ! 215 CALL read_restart() 216 ! 217 ! ... consistency between the input value of nstep and the value 218 ! ... of nstep_path read from the restart_file is checked 219 ! 220 IF ( nstep_path_ == 0 ) THEN 221 ! 222 istep_path = 0 223 nstep_path = nstep_path_ 224 ! 225 END IF 226 ! 227 IF ( nstep_path_ > nstep_path ) nstep_path = nstep_path_ 228 ! 229 ! ... in case first_last_opt has been set to true, reset the frozen 230 ! ... array to false (all the images have to be optimized, at least 231 ! ... on the first iteration) 232 ! 233 IF ( first_last_opt ) frozen = .FALSE. 234 ! 235 ! ... path length is computed here 236 ! 237 path_length = 0.0_DP 238 ! 239 DO i = 1, ( num_of_images - 1 ) 240 ! 241 path_length = path_length + norm( pos(:,i+1) - pos(:,i) ) 242 ! 243 END DO 244 ! 245 ELSE 246 ! 247 CALL initial_guess() 248 ! 249 posold(:,:) = pos(:,:) 250 ! 251 END IF 252 ! 253 ! ... the actual number of degrees of freedom is computed 254 ! 255 deg_of_freedom = 0 256 ! 257 DO i = 1, nat 258 ! 259 IF ( fix_atom_pos(1,i) == 1 ) deg_of_freedom = deg_of_freedom + 1 260 IF ( fix_atom_pos(2,i) == 1 ) deg_of_freedom = deg_of_freedom + 1 261 IF ( fix_atom_pos(3,i) == 1 ) deg_of_freedom = deg_of_freedom + 1 262 ! 263 END DO 264 ! 265 RETURN 266 ! 267 END SUBROUTINE initialize_path 268 ! 269 !-------------------------------------------------------------------- 270 SUBROUTINE initial_guess() 271 !-------------------------------------------------------------------- 272 ! 273 ! ... linear interpolation 274 ! 275 USE path_input_parameters_module, ONLY : input_images 276 USE path_variables, ONLY : pos, dim1, num_of_images, path_length 277 USE path_io_units_module, ONLY : iunpath 278 ! 279 IMPLICIT NONE 280 ! 281 REAL(DP) :: s 282 INTEGER :: i, j 283 LOGICAL :: tooclose 284 REAL(DP), ALLOCATABLE :: pos_n(:,:), dr(:,:), image_spacing(:) 285 ! 286 ! 287 IF ( meta_ionode ) THEN 288 ! 289 ALLOCATE( pos_n( dim1, num_of_images ) ) 290 ALLOCATE( dr( dim1, input_images - 1 ) ) 291 ALLOCATE( image_spacing( input_images - 1 ) ) 292 ! 293 tooclose = .false. 294 image_spacing(:) = 0.0_dp 295 DO i = 1, input_images - 1 296 ! 297 dr(:,i) = ( pos(:,i+1) - pos(:,i) ) 298 ! 299 image_spacing(i) = norm( dr(:,i) ) 300 tooclose = tooclose .OR. ( image_spacing(i) < 0.01 ) 301 ! 302 END DO 303 IF ( tooclose) CALL errore ('initial_guess', & 304 ' something wrong: images are too close',1) 305 ! 306 path_length = SUM( image_spacing(:) ) 307 ! 308 DO i = 1, input_images - 1 309 ! 310 dr(:,i) = dr(:,i) / image_spacing(i) 311 ! 312 END DO 313 ! 314 pos_n(:,1) = pos(:,1) 315 ! 316 i = 1 317 s = 0.0_DP 318 ! 319 DO j = 2, num_of_images - 1 320 ! 321 s = s + path_length / DBLE( num_of_images - 1 ) 322 ! 323 IF ( s > image_spacing(i) ) THEN 324 ! 325 s = s - image_spacing(i) 326 ! 327 i = i + 1 328 ! 329 END IF 330 ! 331 IF ( i >= input_images ) & 332 CALL errore( 'initialize_path', 'i >= input_images', i ) 333 ! 334 pos_n(:,j) = pos(:,i) + s * dr(:,i) 335 ! 336 END DO 337 ! 338 pos_n(:,num_of_images) = pos(:,input_images) 339 ! 340 pos(:,:) = pos_n(:,:) 341 ! 342 path_length = 0.0_DP 343 ! 344 DO i = 1, num_of_images - 1 345 ! 346 path_length = path_length + norm( pos(:,i+1) - pos(:,i) ) 347 ! 348 END DO 349 ! 350 WRITE( UNIT = iunpath, & 351 FMT = '(/,5X,"initial path length",& 352 & T35," = ",F7.4," bohr")' ) path_length 353 ! 354 WRITE( UNIT = iunpath, & 355 FMT = '(5X,"initial inter-image distance",T35," = ",F7.4, & 356 &" bohr")' ) path_length / DBLE( num_of_images - 1 ) 357 ! 358 DEALLOCATE( image_spacing, dr, pos_n ) 359 ! 360 END IF 361 ! 362 CALL mp_bcast( pos, meta_ionode_id, world_comm ) 363 CALL mp_bcast( path_length, meta_ionode_id, world_comm ) 364 ! 365 RETURN 366 ! 367 END SUBROUTINE initial_guess 368 ! 369 !----------------------------------------------------------------------- 370 FUNCTION real_space_tangent( i ) RESULT( rtan ) 371 !----------------------------------------------------------------------- 372 ! 373 ! ... improved definition of the tangent (see JCP 113, 9978) 374 ! 375 USE path_variables, ONLY : dim1, pos, num_of_images, pes 376 ! 377 IMPLICIT NONE 378 ! 379 INTEGER, INTENT(IN) :: i 380 REAL(DP) :: rtan( dim1 ) 381 ! 382 REAL(DP) :: V_previous, V_actual, V_next 383 REAL(DP) :: abs_next, abs_previous 384 REAL(DP) :: delta_V_max, delta_V_min 385 ! 386 ! 387 IF ( i == 1 ) THEN 388 ! 389 rtan(:) = pos(:,i+1) - pos(:,i) 390 ! 391 RETURN 392 ! 393 ELSE IF ( i == num_of_images ) THEN 394 ! 395 rtan(:) = pos(:,i) - pos(:,i-1) 396 ! 397 RETURN 398 ! 399 END IF 400 ! 401 V_previous = pes( i - 1 ) 402 V_actual = pes( i ) 403 V_next = pes( i + 1 ) 404 ! 405 IF ( ( V_next > V_actual ) .AND. ( V_actual > V_previous ) ) THEN 406 ! 407 rtan(:) = pos(:,i+1) - pos(:,i) 408 ! 409 ELSE IF ( ( V_next < V_actual ) .AND. ( V_actual < V_previous ) ) THEN 410 ! 411 rtan(:) = pos(:,i) - pos(:,i-1) 412 ! 413 ELSE 414 ! 415 abs_next = ABS( V_next - V_actual ) 416 abs_previous = ABS( V_previous - V_actual ) 417 ! 418 delta_V_max = MAX( abs_next, abs_previous ) 419 delta_V_min = MIN( abs_next, abs_previous ) 420 ! 421 IF ( V_next > V_previous ) THEN 422 ! 423 rtan(:) = ( pos(:,i+1) - pos(:,i) ) * delta_V_max + & 424 ( pos(:,i) - pos(:,i-1) ) * delta_V_min 425 ! 426 ELSE IF ( V_next < V_previous ) THEN 427 ! 428 rtan(:) = ( pos(:,i+1) - pos(:,i) ) * delta_V_min + & 429 ( pos(:,i) - pos(:,i-1) ) * delta_V_max 430 ! 431 ELSE 432 ! 433 rtan(:) = pos(:,i+1) - pos(:,i-1) 434 ! 435 END IF 436 ! 437 END IF 438 ! 439 rtan(:) = rtan(:) / norm( rtan(:) ) 440 ! 441 RETURN 442 ! 443 END FUNCTION real_space_tangent 444 ! 445 !------------------------------------------------------------------------ 446 SUBROUTINE elastic_constants() 447 !------------------------------------------------------------------------ 448 ! 449 USE path_variables, ONLY : num_of_images, Emax, Emin, & 450 k_max, k_min, k, pes 451 ! 452 IMPLICIT NONE 453 ! 454 INTEGER :: i 455 REAL(DP) :: delta_E 456 REAL(DP) :: k_sum, k_diff 457 ! 458 ! 459 ! ... standard neb ( with springs ) 460 ! 461 k_sum = k_max + k_min 462 k_diff = k_max - k_min 463 ! 464 k(:) = k_min 465 ! 466 delta_E = Emax - Emin 467 ! 468 IF ( delta_E > eps32 ) THEN 469 ! 470 DO i = 1, num_of_images 471 ! 472 k(i) = 0.5_DP*( k_sum - k_diff * & 473 COS( pi * ( pes(i) - Emin ) / delta_E ) ) 474 ! 475 END DO 476 ! 477 END IF 478 ! 479 k(:) = 0.5_DP*k(:) 480 ! 481 RETURN 482 ! 483 END SUBROUTINE elastic_constants 484 ! 485 !------------------------------------------------------------------------ 486 SUBROUTINE neb_gradient() 487 !------------------------------------------------------------------------ 488 ! 489 USE path_variables, ONLY : pos, grad, elastic_grad, grad_pes, k, & 490 num_of_images, climbing, mass, tangent 491 ! 492 IMPLICIT NONE 493 ! 494 INTEGER :: i 495 ! 496 ! 497 IF ( meta_ionode ) THEN 498 ! 499 CALL elastic_constants() 500 ! 501 gradient_loop: DO i = 1, num_of_images 502 ! 503 IF ( i > 1 .AND. i < num_of_images ) THEN 504 ! 505 ! ... elastic gradient only along the path ( variable elastic 506 ! ... consatnt is used ) NEB recipe 507 ! 508 elastic_grad = tangent(:,i) * 0.5_DP * & 509 ( ( k(i) + k(i-1) ) * norm( pos(:,i) - pos(:,(i-1)) ) - & 510 ( k(i) + k(i+1) ) * norm( pos(:,(i+1)) - pos(:,i) ) ) 511 ! 512 END IF 513 ! 514 ! ... total gradient on each image ( climbing image is used if 515 ! ... required ) only the component of the pes gradient orthogonal 516 ! ... to the path is used 517 ! 518 grad(:,i) = grad_pes(:,i) / SQRT( mass(:) ) 519 ! 520 IF ( climbing(i) ) THEN 521 ! 522 grad(:,i) = grad(:,i) - & 523 2.0_DP*tangent(:,i)*( grad(:,i) .dot. tangent(:,i) ) 524 ! 525 ELSE IF ( i > 1 .AND. i < num_of_images ) THEN 526 ! 527 grad(:,i) = elastic_grad + grad(:,i) - & 528 tangent(:,i)*( grad(:,i) .dot. tangent(:,i) ) 529 ! 530 END IF 531 ! 532 END DO gradient_loop 533 ! 534 END IF 535 ! 536 CALL mp_bcast( grad, meta_ionode_id, world_comm ) 537 ! 538 RETURN 539 ! 540 END SUBROUTINE neb_gradient 541 ! 542 !----------------------------------------------------------------------- 543 SUBROUTINE smd_gradient() 544 !----------------------------------------------------------------------- 545 ! 546 USE ions_base, ONLY : if_pos 547 USE path_variables, ONLY : dim1, mass, num_of_images, grad_pes, & 548 tangent, llangevin, lang, grad, ds, & 549 temp_req 550 USE path_variables, ONLY : climbing 551 USE random_numbers, ONLY : gauss_dist 552 ! 553 IMPLICIT NONE 554 ! 555 INTEGER :: i 556 ! 557 ! 558 IF ( meta_ionode ) THEN 559 ! 560 grad(:,:) = 0.0_DP 561 lang(:,:) = 0.0_DP 562 ! 563 ! ... we project pes gradients and gaussian noise 564 ! 565 DO i = 1, num_of_images 566 ! 567 IF ( llangevin ) THEN 568 ! 569 ! ... the random term used in langevin dynamics is generated here 570 ! 571 lang(:,i) = gauss_dist( 0.0_DP, SQRT( 2.0_DP*temp_req*ds ), dim1 ) 572 ! 573 lang(:,i) = lang(:,i)*DBLE( RESHAPE( if_pos, (/ dim1 /) ) ) 574 ! 575 END IF 576 ! 577 grad(:,i) = grad_pes(:,i) / SQRT( mass(:) ) 578 ! 579 IF ( climbing(i) ) THEN 580 ! 581 grad(:,i) = grad(:,i) - & 582 2.0_DP*tangent(:,i)*( grad(:,i) .dot. tangent(:,i) ) 583 ! 584 ELSE IF ( i > 1 .AND. i < num_of_images ) THEN 585 ! 586 ! ... projection of the pes gradients 587 ! 588 grad(:,i) = grad(:,i) - & 589 tangent(:,i)*( grad(:,i) .dot. tangent(:,i) ) 590 ! 591 IF ( llangevin ) THEN 592 ! 593 lang(:,i) = lang(:,i) - & 594 tangent(:,i)*( lang(:,i) .dot. tangent(:,i) ) 595 ! 596 END IF 597 ! 598 END IF 599 ! 600 END DO 601 ! 602 END IF 603 ! 604 CALL mp_bcast( grad, meta_ionode_id, world_comm ) 605 CALL mp_bcast( lang, meta_ionode_id, world_comm ) 606 ! 607 RETURN 608 ! 609 END SUBROUTINE smd_gradient 610 ! 611 ! ... shared routines 612 ! 613 !----------------------------------------------------------------------- 614 FUNCTION new_tangent() RESULT( ntan ) 615 !----------------------------------------------------------------------- 616 ! 617 USE path_variables, ONLY : dim1, num_of_images 618 ! 619 IMPLICIT NONE 620 ! 621 REAL(DP) :: ntan( dim1, num_of_images ) 622 ! 623 INTEGER :: i 624 ! 625 ! 626 IF ( meta_ionode ) THEN 627 ! 628 DO i = 1, num_of_images 629 ! 630 ntan(:,i) = real_space_tangent( i ) 631 ! 632 END DO 633 ! 634 END IF 635 ! 636 CALL mp_bcast( ntan, meta_ionode_id, world_comm ) 637 ! 638 RETURN 639 ! 640 END FUNCTION new_tangent 641 ! 642 !----------------------------------------------------------------------- 643 SUBROUTINE compute_error( err_out ) 644 !----------------------------------------------------------------------- 645 ! 646 USE path_variables, ONLY : pos, posold, num_of_images, grad, & 647 use_freezing, first_last_opt, path_thr, & 648 error, frozen, lquick_min 649 USE mp_images, ONLY : nimage 650 ! 651 IMPLICIT NONE 652 ! 653 REAL(DP), OPTIONAL, INTENT(OUT) :: err_out 654 ! 655 INTEGER :: i 656 INTEGER :: fii, lii, freed, num_of_scf_images 657 REAL(DP) :: err_max 658 LOGICAL :: first 659 ! 660 ! 661 IF ( first_last_opt ) THEN 662 ! 663 fii = 1 664 lii = num_of_images 665 ! 666 frozen = .FALSE. 667 ! 668 ELSE 669 ! 670 fii = 2 671 lii = num_of_images - 1 672 ! 673 frozen = .FALSE. 674 ! 675 ! ... the first and the last images are always frozen 676 ! 677 frozen(fii-1) = .TRUE. 678 frozen(lii+1) = .TRUE. 679 ! 680 END IF 681 ! 682 IF ( meta_ionode ) THEN 683 ! 684 DO i = 1, num_of_images 685 ! 686 ! ... the error is given by the largest component of the gradient 687 ! ... vector ( PES + SPRINGS in the neb case ) 688 ! 689 error(i) = MAXVAL( ABS( grad(:,i) ) ) / bohr_radius_angs * autoev 690 ! 691 END DO 692 ! 693 err_max = MAXVAL( error(fii:lii), 1 ) 694 ! 695 IF ( use_freezing ) THEN 696 ! 697 frozen(fii:lii) = ( error(fii:lii) < & 698 MAX( 0.5_DP*err_max, path_thr ) ) 699 ! 700 END IF 701 ! 702 IF ( nimage > 1 .AND. use_freezing ) THEN 703 ! 704 find_scf_images: DO 705 ! 706 num_of_scf_images = COUNT( .NOT.frozen(fii:lii) ) 707 ! 708 IF ( num_of_scf_images >= nimage ) EXIT find_scf_images 709 ! 710 first = .TRUE. 711 ! 712 search: DO i = fii, lii 713 ! 714 IF ( .NOT.frozen(i) ) CYCLE search 715 ! 716 IF ( first ) THEN 717 ! 718 first = .FALSE. 719 freed = i 720 ! 721 CYCLE search 722 ! 723 END IF 724 ! 725 IF ( error(i) > error(freed) ) freed = i 726 ! 727 END DO search 728 ! 729 frozen(freed) = .FALSE. 730 ! 731 END DO find_scf_images 732 ! 733 END IF 734 ! 735 IF ( use_freezing .AND. lquick_min ) THEN 736 ! 737 ! ... the old positions of the frozen images are set to the 738 ! ... present position (equivalent to resetting the velocity) 739 ! 740 FORALL( i = fii:lii, frozen(i) ) posold(:,i) = pos(:,i) 741 ! 742 END IF 743 ! 744 END IF 745 ! 746 CALL mp_bcast( error, meta_ionode_id, world_comm ) 747 CALL mp_bcast( err_max, meta_ionode_id, world_comm ) 748 CALL mp_bcast( frozen, meta_ionode_id, world_comm ) 749 CALL mp_bcast( posold, meta_ionode_id, world_comm ) 750 ! 751 IF ( PRESENT( err_out ) ) err_out = err_max 752 ! 753 RETURN 754 ! 755 END SUBROUTINE compute_error 756 ! 757 !----------------------------------------------------------------------- 758 SUBROUTINE fcp_compute_error( err_out ) 759 !----------------------------------------------------------------------- 760 ! 761 USE path_variables, ONLY : num_of_images, first_last_opt 762 USE fcp_variables, ONLY : fcp_mu 763 USE fcp_opt_routines, ONLY : fcp_neb_ef 764 ! 765 IMPLICIT NONE 766 ! 767 REAL(DP), OPTIONAL, INTENT(OUT) :: err_out 768 ! 769 INTEGER :: i 770 INTEGER :: fii, lii 771 REAL(DP) :: err_max 772 ! 773 ! 774 IF ( first_last_opt ) THEN 775 ! 776 fii = 1 777 lii = num_of_images 778 ! 779 ELSE 780 ! 781 fii = 2 782 lii = num_of_images - 1 783 ! 784 END IF 785 ! 786 IF ( meta_ionode ) THEN 787 ! 788 err_max = MAXVAL( ABS( fcp_mu - fcp_neb_ef(fii:lii) ), 1 ) 789 ! 790 END IF 791 ! 792 CALL mp_bcast( err_max, meta_ionode_id, world_comm ) 793 ! 794 IF ( PRESENT( err_out ) ) err_out = err_max 795 ! 796 RETURN 797 ! 798 END SUBROUTINE fcp_compute_error 799 ! 800 !------------------------------------------------------------------------ 801 SUBROUTINE born_oppenheimer_pes( stat ) 802 !------------------------------------------------------------------------ 803 ! 804 USE path_variables, ONLY : num_of_images, & 805 pending_image, istep_path, pes, & 806 first_last_opt, Emin, Emax, Emax_index 807 ! 808 IMPLICIT NONE 809 ! 810 LOGICAL, INTENT(OUT) :: stat 811 ! 812 INTEGER :: fii, lii 813 ! 814 ! 815 IF ( istep_path == 0 .OR. first_last_opt ) THEN 816 ! 817 fii = 1 818 lii = num_of_images 819 ! 820 ELSE 821 ! 822 fii = 2 823 lii = num_of_images - 1 824 ! 825 END IF 826 ! 827 IF ( pending_image /= 0 ) fii = pending_image 828 ! 829 CALL compute_scf( fii, lii, stat ) 830 ! 831 IF ( .NOT. stat ) RETURN 832 ! 833 Emin = MINVAL( pes(1:num_of_images) ) 834 Emax = MAXVAL( pes(1:num_of_images) ) 835 Emax_index = MAXLOC( pes(1:num_of_images), 1 ) 836 ! 837 RETURN 838 ! 839 END SUBROUTINE born_oppenheimer_pes 840 ! 841 !------------------------------------------------------------------------ 842 SUBROUTINE fe_profile() 843 !------------------------------------------------------------------------ 844 ! 845 USE path_variables, ONLY : num_of_images 846 USE path_variables, ONLY : pos, pes, grad_pes, & 847 Emin, Emax, Emax_index 848 ! 849 IMPLICIT NONE 850 ! 851 INTEGER :: i 852 ! 853 ! 854 pes(:) = 0.0_DP 855 ! 856 DO i = 2, num_of_images 857 ! 858 pes(i) = pes(i-1) + 0.5_DP*( ( pos(:,i) - pos(:,i-1) ) .dot. & 859 ( grad_pes(:,i) + grad_pes(:,i-1) ) ) 860 ! 861 END DO 862 ! 863 Emin = MINVAL( pes(1:num_of_images) ) 864 Emax = MAXVAL( pes(1:num_of_images) ) 865 Emax_index = MAXLOC( pes(1:num_of_images), 1 ) 866 ! 867 RETURN 868 ! 869 END SUBROUTINE fe_profile 870 ! 871 !----------------------------------------------------------------------- 872 SUBROUTINE search_mep() 873 !----------------------------------------------------------------------- 874 ! 875 USE path_variables, ONLY : lneb, lsmd 876 USE path_variables, ONLY : conv_path, istep_path, nstep_path, & 877 pending_image, activation_energy, & 878 err_max, pes, climbing, CI_scheme, & 879 Emax_index, fixed_tan, tangent 880 USE path_io_routines, ONLY : write_restart, write_dat_files, write_output 881 USE path_formats, ONLY : scf_iter_fmt 882 USE fcp_variables, ONLY : lfcpopt 883 ! 884 USE path_reparametrisation 885 ! 886 IMPLICIT NONE 887 ! 888 LOGICAL :: stat 889 REAL(DP) :: fcp_err_max = 0.0_DP 890 ! 891 REAL(DP), EXTERNAL :: get_clock 892 ! 893 ! 894 conv_path = .FALSE. 895 ! 896 CALL search_mep_init() 897 ! 898 IF ( istep_path == nstep_path ) THEN 899 ! 900 CALL write_dat_files() 901 ! 902 CALL write_output() 903 ! 904 pending_image = 0 905 ! 906 CALL write_restart() 907 ! 908 RETURN 909 ! 910 END IF 911 ! 912 ! ... path optimisation loop 913 ! 914 optimisation: DO 915 ! 916 ! ... new positions are saved on file: it has to be done here 917 ! ... because, in the event of an unexpected crash the new positions 918 ! ... would be lost. At this stage the forces and the energies are 919 ! ... not yet known (but are not necessary for restarting); the 920 ! ... restart file is written again as soon as the energies and 921 ! ... forces have been computed. 922 ! 923 CALL write_restart() 924 ! 925 IF ( meta_ionode ) & 926 WRITE( UNIT = iunpath, FMT = scf_iter_fmt ) istep_path + 1 927 ! 928 ! ... energies and gradients acting on each image of the path (in real 929 ! ... space) are computed calling a driver for the scf calculations 930 ! 931 ! 932 CALL born_oppenheimer_pes( stat ) 933 ! 934 ! 935 IF ( .NOT. stat ) THEN 936 ! 937 conv_path = .FALSE. 938 ! 939 EXIT optimisation 940 ! 941 END IF 942 ! 943 ! ... istep_path is updated after a self-consistency step has been 944 ! ... completed 945 ! 946 istep_path = istep_path + 1 947 ! 948 ! ... the new tangent is computed here : 949 ! ... the improved definition of tangent requires the energies 950 ! 951 IF ( .NOT. fixed_tan ) tangent(:,:) = new_tangent() 952 ! 953 IF ( CI_scheme == "auto" ) THEN 954 ! 955 climbing = .FALSE. 956 ! 957 climbing(Emax_index) = .TRUE. 958 ! 959 END IF 960 ! 961 IF ( lneb ) CALL neb_gradient() 962 IF ( lsmd ) CALL smd_gradient() 963 ! 964 ! ... the forward activation energy is computed here 965 ! 966 activation_energy = ( pes(Emax_index) - pes(1) )*autoev 967 ! 968 ! ... the error is computed here (frozen images are also set here) 969 ! 970 CALL compute_error( err_max ) 971 IF ( lfcpopt ) CALL fcp_compute_error( fcp_err_max ) 972 ! 973 ! ... information is written on the files 974 ! 975 CALL write_dat_files() 976 ! 977 ! ... information is written on the standard output 978 ! 979 CALL write_output() 980 ! 981 ! ... the restart file is written 982 ! 983 CALL write_restart() 984 ! 985 ! ... exit conditions 986 ! 987 IF ( check_exit( err_max, fcp_err_max ) ) EXIT optimisation 988 ! 989 ! ... if convergence is not yet achieved, the path is optimised 990 ! 991 CALL optimisation_step() 992 ! 993 IF ( lsmd ) CALL reparametrise() 994 ! 995 END DO optimisation 996 ! 997 ! ... the restart file is written before the exit 998 ! 999 CALL write_restart() 1000 ! 1001 RETURN 1002 ! 1003 END SUBROUTINE search_mep 1004 ! 1005 !------------------------------------------------------------------------ 1006 SUBROUTINE search_mep_init() 1007 !------------------------------------------------------------------------ 1008 ! 1009 USE path_variables, ONLY : lsmd 1010 USE path_variables, ONLY : pending_image, tangent 1011 ! 1012 USE path_reparametrisation 1013 ! 1014 IMPLICIT NONE 1015 ! 1016 ! 1017 IF ( pending_image /= 0 ) RETURN 1018 ! 1019 IF ( lsmd ) CALL reparametrise() 1020 ! 1021 tangent(:,:) = new_tangent() 1022 ! 1023 RETURN 1024 ! 1025 END SUBROUTINE search_mep_init 1026 ! 1027 !------------------------------------------------------------------------ 1028 FUNCTION check_exit( err_max, fcp_err_max ) 1029 !------------------------------------------------------------------------ 1030 ! 1031 USE path_input_parameters_module, ONLY : num_of_images_inp => num_of_images 1032 USE path_variables, ONLY : lneb, lsmd 1033 USE path_variables, ONLY : path_thr, istep_path, nstep_path, & 1034 conv_path, pending_image, & 1035 num_of_images, llangevin 1036 USE path_formats, ONLY : final_fmt 1037 USE fcp_variables, ONLY : lfcpopt, fcp_relax_crit 1038 ! 1039 IMPLICIT NONE 1040 ! 1041 LOGICAL :: check_exit 1042 REAL(DP), INTENT(IN) :: err_max 1043 REAL(DP), INTENT(IN) :: fcp_err_max 1044 LOGICAL :: exit_condition 1045 ! 1046 ! 1047 check_exit = .FALSE. 1048 ! 1049 ! ... the program checks if the convergence has been achieved 1050 ! 1051 exit_condition = ( .NOT.llangevin .AND. & 1052 ( num_of_images == num_of_images_inp ) .AND. & 1053 ( err_max <= path_thr ) ) 1054 ! 1055 IF ( lfcpopt .AND. fcp_err_max > fcp_relax_crit ) THEN 1056 exit_condition = .FALSE. 1057 END IF 1058 ! 1059 IF ( exit_condition ) THEN 1060 ! 1061 IF ( meta_ionode ) THEN 1062 ! 1063 WRITE( UNIT = iunpath, FMT = final_fmt ) 1064 ! 1065 IF ( lneb ) & 1066 WRITE( UNIT = iunpath, & 1067 FMT = '(/,5X,"neb: convergence achieved in ",I3, & 1068 & " iterations" )' ) istep_path 1069 IF ( lsmd ) & 1070 WRITE( UNIT = iunpath, & 1071 FMT = '(/,5X,"smd: convergence achieved in ",I3, & 1072 & " iterations" )' ) istep_path 1073 ! 1074 END IF 1075 ! 1076 pending_image = 0 1077 ! 1078 conv_path = .TRUE. 1079 check_exit = .TRUE. 1080 ! 1081 RETURN 1082 ! 1083 END IF 1084 ! 1085 ! ... the program checks if the maximum number of iterations has 1086 ! ... been reached 1087 ! 1088 IF ( istep_path >= nstep_path ) THEN 1089 ! 1090 IF ( meta_ionode ) THEN 1091 ! 1092 WRITE( UNIT = iunpath, FMT = final_fmt ) 1093 ! 1094 IF ( lneb ) & 1095 WRITE( UNIT = iunpath, & 1096 FMT = '(/,5X,"neb: reached the maximum number of ", & 1097 & "steps")' ) 1098 IF ( lsmd ) & 1099 WRITE( UNIT = iunpath, & 1100 FMT = '(/,5X,"smd: reached the maximum number of ", & 1101 & "steps")' ) 1102 ! 1103 END IF 1104 ! 1105 pending_image = 0 1106 ! 1107 check_exit = .TRUE. 1108 ! 1109 RETURN 1110 ! 1111 END IF 1112 ! 1113 RETURN 1114 ! 1115 END FUNCTION check_exit 1116 ! 1117 !------------------------------------------------------------------------ 1118 SUBROUTINE optimisation_step() 1119 !------------------------------------------------------------------------ 1120 ! 1121 USE path_variables, ONLY : num_of_images, frozen, lsteep_des, & 1122 lquick_min, lbroyden, lbroyden2, & 1123 llangevin, istep_path 1124 USE path_opt_routines, ONLY : quick_min, broyden, broyden2, & 1125 steepest_descent, langevin 1126 USE fcp_variables, ONLY : lfcpopt 1127 USE fcp_opt_routines, ONLY : fcp_opt_perform 1128 ! 1129 IMPLICIT NONE 1130 ! 1131 INTEGER :: image 1132 ! 1133 ! 1134 IF ( lbroyden ) THEN 1135 ! 1136 CALL broyden() 1137 ! 1138 ELSE IF (lbroyden2 ) THEN 1139 ! 1140 CALL broyden2() 1141 ! 1142 ELSE 1143 ! 1144 DO image = 1, num_of_images 1145 ! 1146 IF ( frozen(image) ) CYCLE 1147 ! 1148 IF ( lsteep_des ) THEN 1149 ! 1150 CALL steepest_descent( image ) 1151 ! 1152 ELSE IF ( llangevin ) THEN 1153 ! 1154 CALL langevin( image ) 1155 ! 1156 ELSE IF ( lquick_min ) THEN 1157 ! 1158 CALL quick_min( image, istep_path ) 1159 ! 1160 END IF 1161 ! 1162 END DO 1163 ! 1164 END IF 1165 ! 1166 IF ( lfcpopt ) CALL fcp_opt_perform() 1167 ! 1168 RETURN 1169 ! 1170 END SUBROUTINE optimisation_step 1171 ! 1172END MODULE path_base 1173