1! 2! Copyright (C) 2002-2010 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 SUBROUTINE ecutoffs_setup( ecutwfc_, ecutrho_, ecfixed_, qcutz_, & 10 q2sigma_, refg_ ) 11!------------------------------------------------------------------------------! 12 USE kinds, ONLY: DP 13 USE constants, ONLY: eps8 14 USE gvecw, ONLY: ecutwfc 15 USE gvecw, ONLY: ecfixed, qcutz, q2sigma 16 USE gvect, ONLY: ecutrho 17 USE gvecs, ONLY: ecuts, dual, doublegrid 18 USE pseudopotential, only: tpstab 19 USE io_global, only: stdout, ionode 20 USE uspp, only: okvan 21 use betax, only: mmx, refg 22 23 IMPLICIT NONE 24 REAL(DP), INTENT(IN) :: ecutwfc_, ecutrho_, ecfixed_, qcutz_, & 25 q2sigma_, refg_ 26 27 ecutwfc = ecutwfc_ 28 29 IF ( ecutrho_ <= 0.D0 ) THEN 30 ! 31 dual = 4.D0 32 ! 33 ELSE 34 ! 35 dual = ecutrho_ / ecutwfc 36 ! 37 IF ( dual <= 1.D0 ) & 38 CALL errore( ' ecutoffs_setup ', ' invalid dual? ', 1 ) 39 ! 40 END IF 41 42 doublegrid = ( dual > 4.0_dp + eps8 ) 43 IF ( doublegrid .AND. .NOT. okvan ) & 44 CALL errore( 'setup', 'No USPP: set ecutrho=4*ecutwfc', 1 ) 45 ecutrho = dual * ecutwfc 46 ! 47 IF ( doublegrid ) THEN 48 ! 49 ecuts = 4.D0 * ecutwfc 50 ! 51 ELSE 52 ! 53 ecuts = ecutrho 54 ! 55 END IF 56 ! 57 ecfixed = ecfixed_ 58 qcutz = qcutz_ 59 q2sigma = q2sigma_ 60 61 IF( refg_ < 0.0001d0 ) THEN 62 tpstab = .FALSE. 63 refg = 0.05d0 64 ELSE 65 refg = refg_ 66 END IF 67 68 CALL set_interpolation_table_size( mmx, refg, ecutrho ) 69 70 RETURN 71 END SUBROUTINE ecutoffs_setup 72 73 74 SUBROUTINE set_interpolation_table_size( mmx, refg, gmax ) 75 USE control_flags, only: thdyn 76 USE kinds, only: DP 77 IMPLICIT NONE 78 INTEGER, INTENT(OUT) :: mmx 79 REAL(DP), INTENT(IN) :: refg 80 REAL(DP), INTENT(IN) :: gmax 81 IF( thdyn ) THEN 82 ! ... a larger table is used when cell is moving to allow 83 ! ... large volume fluctuation 84 mmx = NINT( 2.0d0 * gmax / refg ) 85 ELSE 86 mmx = NINT( 1.2d0 * gmax / refg ) 87 END IF 88 RETURN 89 END SUBROUTINE set_interpolation_table_size 90 91 92 SUBROUTINE gcutoffs_setup( alat, tk_inp, nk_inp, kpoints_inp ) 93 94! (describe briefly what this routine does...) 95! ---------------------------------------------- 96 97 USE kinds, ONLY: DP 98 USE gvecw, ONLY: ecutwfc, gcutw 99 USE gvect, ONLY: ecutrho, gcutm 100 USE gvecs, ONLY: ecuts, gcutms 101 USE gvecw, ONLY: ekcut, gkcut 102 USE constants, ONLY: eps8, pi 103 104 IMPLICIT NONE 105 106! ... declare subroutine arguments 107 REAL(DP), INTENT(IN) :: alat 108 LOGICAL, INTENT(IN) :: tk_inp 109 INTEGER, INTENT(IN) :: nk_inp 110 REAL(DP), INTENT(IN) :: kpoints_inp(3,*) 111 112! ... declare other variables 113 INTEGER :: i 114 REAL(DP) :: kcut, ksq 115 REAL(DP) :: tpiba 116 117! end of declarations 118! ---------------------------------------------- 119 120! ... Set Values for the cutoff 121 122 123 IF( alat < eps8 ) THEN 124 CALL errore(' cut-off setup ', ' alat too small ', 0) 125 END IF 126 127 tpiba = 2.0d0 * pi / alat 128 129 ! ... Constant cutoff simulation parameters 130 131 gcutw = ecutwfc / tpiba**2 ! wave function cut-off 132 gcutm = ecutrho / tpiba**2 ! potential cut-off 133 gcutms= ecuts / tpiba**2 ! smooth mesh cut-off 134 135 kcut = 0.0_DP 136 IF ( tk_inp ) THEN 137! ... augment plane wave cutoff to include all k+G's 138 DO i = 1, nk_inp 139! ... calculate modulus 140 ksq = kpoints_inp( 1, i ) ** 2 + kpoints_inp( 2, i ) ** 2 + kpoints_inp( 3, i ) ** 2 141 IF ( ksq > kcut ) kcut = ksq 142 END DO 143 END IF 144 145 gkcut = ( sqrt( kcut ) + sqrt( gcutw ) ) ** 2 146 147 ekcut = gkcut * tpiba ** 2 148 149 RETURN 150 END SUBROUTINE gcutoffs_setup 151 152! ---------------------------------------------- 153 154 SUBROUTINE cutoffs_print_info() 155 156 ! Print out information about different cut-offs 157 158 USE gvecw, ONLY: ecutwfc, gcutw 159 USE gvect, ONLY: ecutrho, gcutm 160 USE gvecw, ONLY: ecfixed, qcutz, q2sigma 161 USE gvecw, ONLY: ekcut, gkcut 162 USE gvecs, ONLY: ecuts, gcutms 163 use betax, only: mmx, refg 164 USE io_global, ONLY: stdout 165 USE input_parameters, ONLY: ref_cell, ref_alat 166 167 WRITE( stdout, 100 ) ecutwfc, ecutrho, ecuts, sqrt(gcutw), & 168 sqrt(gcutm), sqrt(gcutms) 169 170 IF(ref_cell) WRITE( stdout,'(3X,"Reference Cell alat is",F14.8,1X,"A.U is used to Compute Gcutoffs:")') ref_alat ! BS : debug 171 172 IF( qcutz > 0.0d0 ) THEN 173 WRITE( stdout, 150 ) qcutz, q2sigma, ecfixed 174 END IF 175 176 WRITE( stdout,200) refg, mmx 177 178100 FORMAT(/,3X,'Energy Cut-offs',/ & 179 ,3X,'---------------',/ & 180 ,3X,'Ecutwfc = ',F6.1,' Ry, ', 3X,'Ecutrho = ',F6.1,' Ry, ', 3X,'Ecuts = ',F6.1,' Ry',/ & 181 ,3X,'Gcutwfc = ',F6.1,' , ', 3X,'Gcutrho = ',F6.1,' ', 3X,'Gcuts = ',F6.1) 182150 FORMAT( 3X,'modified kinetic energy functional, with parameters:',/, & 183 3X,'ecutz = ',f8.4,' ecsig = ', f7.4,' ecfix = ',f6.2) 184200 FORMAT( 3X,'NOTA BENE: refg, mmx = ', f10.6,I6 ) 185 186 RETURN 187 END SUBROUTINE cutoffs_print_info 188 189! ---------------------------------------------- 190 191 SUBROUTINE orthogonalize_info( ) 192 USE control_flags, ONLY: ortho_eps, ortho_max 193 USE io_global, ONLY: stdout 194 IMPLICIT NONE 195 WRITE(stdout, 585) 196 WRITE(stdout, 511) ortho_eps, ortho_max 197 511 FORMAT( 3X,'Orthog. with lagrange multipliers : eps = ',E10.2, ', max = ',I3) 198 585 FORMAT( 3X,'Eigenvalues calculated without the kinetic term contribution') 199 RETURN 200 END SUBROUTINE orthogonalize_info 201 202 203! ---------------------------------------------- 204 205 206 SUBROUTINE electrons_print_info( ) 207 208 USE kinds, ONLY: DP 209 USE electrons_base, ONLY: nbnd, nspin, nel, nelt, nupdwn, iupdwn, & 210 f, qbac 211 USE io_global, ONLY: stdout 212 USE ions_base, ONLY: zv, nsp, na 213 214 IMPLICIT NONE 215 INTEGER :: i,is 216 217 IF( nspin == 1) THEN 218 WRITE(stdout,6) nelt, nbnd 219 WRITE(stdout,7) ( f( i ), i = 1, nbnd ) 220 ELSE 221 WRITE(stdout,8) nelt 222 WRITE(stdout,9) nel(1) 223 WRITE(stdout,7) ( f( i ), i = 1, nupdwn(1)) 224 WRITE(stdout,10) nel(2) 225 WRITE(stdout,7) ( f( i ), i = iupdwn(2), ( iupdwn(2) + nupdwn(2) - 1 ) ) 226 END IF 227 228 qbac=0. 229 do is=1,nsp 230 qbac=qbac+na(is)*zv(is) 231 end do 232 qbac=qbac-nelt 233 if(qbac.ne.0) write(stdout,11) qbac 234 235 2366 FORMAT(/,3X,'Electronic states',/ & 237 ,3X,'-----------------',/ & 238 ,3X,'Number of Electrons= ',I5,', of States = ',I5,/ & 239 ,3X,'Occupation numbers :') 2407 FORMAT(2X,10F5.2) 2418 FORMAT(/,3X,'Electronic states',/ & 242 ,3X,'-----------------',/ & 243 ,3X,'Local Spin Density calculation',/ & 244 ,3X,'Number of Electrons= ',I5) 2459 FORMAT( 3X,'Spins up = ', I5, ', occupations: ') 24610 FORMAT( 3X,'Spins down = ', I5, ', occupations: ') 24711 FORMAT(/,3X,'WARNING: system charge = ',F12.6) 248 RETURN 249 END SUBROUTINE electrons_print_info 250 251 252! ---------------------------------------------- 253 254 255 SUBROUTINE exch_corr_print_info() 256 257 USE funct, ONLY: write_dft_name 258 USE io_global, ONLY: stdout 259 260 IMPLICIT NONE 261 262 WRITE(stdout,800) 263 call write_dft_name ( ) 264800 FORMAT(//,3X,'Exchange and correlations functionals',/ & 265 ,3X,'-------------------------------------') 266 267 RETURN 268 END SUBROUTINE exch_corr_print_info 269 270 271 272! ---------------------------------------------- 273 274 275 276 SUBROUTINE ions_print_info( ) 277 278 ! Print info about input parameter for ion dynamic 279 280 USE io_global, ONLY: ionode, stdout 281 USE control_flags, ONLY: tranp, amprp, tnosep, tolp, tfor, tsdp, & 282 tzerop, tv0rd, taurdr, nbeg, tcp, tcap 283 USE ions_base, ONLY: if_pos, nsp, na, tau, ityp, & 284 amass, nat, fricp, greasp, rcmax 285 USE ions_nose, ONLY: tempw, ndega 286 USE constants, ONLY: amu_au 287 288 IMPLICIT NONE 289 290 integer is, ia, k, ic 291 LOGICAL :: ismb( 3 ) 292 293 WRITE( stdout, 50 ) 294 295 IF( .NOT. tfor ) THEN 296 WRITE( stdout, 518 ) 297 ELSE 298 WRITE( stdout, 520 ) 299 IF( tsdp ) THEN 300 WRITE( stdout, 521 ) 301 ELSE 302 WRITE( stdout, 522 ) 303 END IF 304 WRITE( stdout, 523 ) ndega 305 WRITE( stdout, 524 ) fricp, greasp 306 IF( tv0rd ) THEN 307 WRITE( stdout, 850 ) 308 ELSE IF ( tzerop ) THEN 309 WRITE( stdout, 635 ) 310 ENDIF 311 END IF 312 313 DO is = 1, nsp 314 IF( tranp(is) ) THEN 315 WRITE( stdout,510) 316 WRITE( stdout,512) is, amprp(is) 317 END IF 318 END DO 319 320 WRITE(stdout,660) 321 DO is = 1, nsp 322 WRITE(stdout,1000) is, na(is), amass(is)*amu_au, amass(is), rcmax(is) 323 DO ia = 1, nat 324 IF( ityp(ia) == is ) THEN 325 WRITE(stdout,1010) ( tau(k,ia), K = 1,3 ) 326 END IF 327 END DO 328 END DO 329 330 IF ( ( nbeg > -1 ) .AND. ( .NOT. taurdr ) ) THEN 331 WRITE(stdout,661) 332 ELSE 333 WRITE(stdout,662) 334 ENDIF 335 336 IF( tfor ) THEN 337 338 IF( ANY( ( if_pos( 1:3, 1:nat ) == 0 ) ) ) THEN 339 340 WRITE(stdout,1020) 341 WRITE(stdout,1022) 342 343 DO ia = 1, nat 344 ismb( 1 ) = ( if_pos(1,ia) /= 0 ) 345 ismb( 2 ) = ( if_pos(2,ia) /= 0 ) 346 ismb( 3 ) = ( if_pos(3,ia) /= 0 ) 347 IF( .NOT. ALL( ismb ) ) THEN 348 WRITE( stdout, 1023 ) ia, ( ismb(k), K = 1, 3 ) 349 END IF 350 END DO 351 352 ELSE 353 354 WRITE(stdout,1021) 355 356 END IF 357 END IF 358 359 IF( tfor ) THEN 360 if( ( tcp .or. tcap .or. tnosep ) .and. tsdp ) then 361 call errore(' ions_print_info', & 362 ' Temperature control not allowed with steepest descent',1) 363 endif 364 IF(.not. tcp .and. .not. tcap .and. .not. tnosep ) THEN 365 WRITE( stdout,550) 366 ELSE IF( tcp .and. tcap ) then 367 call errore(' ions_print_info', ' Velocity rescaling not' & 368 //' compatible with random velocity initialization',1) 369 ELSE IF( tcp .and. tnosep ) then 370 call errore(' ions_print_info', ' Velocity rescaling and' & 371 //' Nose thermostat are incompatible',1) 372 ELSE IF(tcap .and. tnosep ) then 373 call errore(' ions_print_info', ' Nose thermostat not' & 374 //' compatible with random velocity initialization',1) 375 ELSE IF(tcp) THEN 376 WRITE( stdout,555) tempw,tolp 377 !ELSE IF(tcap) THEN !tcap is random velocity initialization! 378 ! WRITE( stdout,560) tempw,tolp 379 ELSE IF(tnosep) THEN 380 WRITE( stdout,595) 381 ELSE 382 WRITE( stdout,550) 383 END IF 384 END IF 385 386 50 FORMAT(//,3X,'Ions Simulation Parameters',/ & 387 ,3X,'--------------------------') 388 389 510 FORMAT( 3X,'Initial random displacement of ionic coordinates',/, & 390 3X,' specie amplitude') 391 512 FORMAT( 3X,I7,2X,F9.6) 392 393 518 FORMAT( 3X,'Ions are not allowed to move') 394 520 FORMAT( 3X,'Ions are allowed to move') 395 521 FORMAT( 3X,'Ions dynamics with steepest descent') 396 522 FORMAT( 3X,'Ions dynamics with newton equations') 397 523 format( 3X,'the temperature is computed for ',i5,' degrees of freedom') 398 524 format( 3X,'ion dynamics with fricp = ',f7.4,' and greasp = ',f7.4) 399 550 FORMAT( 3X,'Ionic temperature is not controlled') 400 555 FORMAT( 3X,'Ionic temperature control via ', & 401 'rescaling of velocities :',/ & 402 ,3X,'temperature required = ',F10.5,'K, ', & 403 'tolerance = ',F10.5,'K') 404 560 FORMAT( 3X,'Ionic temperature control via ', & 405 'canonical velocities rescaling :',/ & 406 ,3X,'temperature required = ',F10.5,'K, ', & 407 'tolerance = ',F10.5,'K') 408 595 FORMAT( 3X,'Ionic temperature control via nose thermostat') 409 635 FORMAT( 3X,'Zero initial momentum for ions') 410 411 660 FORMAT( 3X,'Ionic position (from input)', /, & 412 3X,'sorted by specie, and converted to real a.u. coordinates') 413 661 FORMAT( 3X,'Ionic position will be re-read from restart file') 414 662 FORMAT( 3X,'Ionic position read from input file') 415 416 850 FORMAT( 3X,'Initial ion velocities read from input') 417 418 1000 FORMAT(3X,'Species ',I3,' atoms = ',I4,' mass = ',F12.2, ' (a.u.), ', & 419 & F12.2, ' (amu)', ' rcmax = ', F6.2, ' (a.u.)' ) 420 1010 FORMAT(3X,3(1X,F12.6)) 421 1020 FORMAT(/,3X,'NOT all atoms are allowed to move ') 422 1021 FORMAT(/,3X,'All atoms are allowed to move') 423 1022 FORMAT( 3X,' indx ..x.. ..y.. ..z..') 424 1023 FORMAT( 3X,I4,3(1X,L5)) 425 426 427 428 RETURN 429 END SUBROUTINE ions_print_info 430 431 432! ---------------------------------------------- 433 434 subroutine cell_print_info( ) 435 436 USE constants, ONLY: au_gpa 437 USE control_flags, ONLY: thdyn, tsdc, tzeroc, tbeg, nbeg, tpre 438 USE control_flags, ONLY: tnoseh 439 USE io_global, ONLY: stdout 440 USE cell_base, ONLY: press, frich, greash, wmass 441 442 IMPLICIT NONE 443 444 WRITE(stdout,545 ) 445 IF ( tpre ) WRITE( stdout, 600 ) 446 IF ( tbeg ) THEN 447 WRITE(stdout,546) 448 ELSE 449 WRITE(stdout,547) 450 IF( nbeg > -1 ) WRITE( stdout, 548 ) 451 END IF 452 453 IF( .NOT. thdyn ) THEN 454 WRITE( stdout,525) 455 WRITE( stdout,606) 456 ELSE 457 IF( tsdc ) THEN 458 WRITE( stdout,526) 459 ELSE 460 IF( frich /= 0.0d0 ) THEN 461 WRITE( stdout,602) frich, greash 462 ELSE 463 WRITE( stdout,527) 464 END IF 465 IF( tnoseh ) then 466 WRITE( stdout,604) 467 ELSE 468 WRITE( stdout,565) 469 END IF 470 IF( tzeroc ) THEN 471 WRITE( stdout,563) 472 ENDIF 473 END IF 474 WRITE( stdout,530) press * au_gpa, wmass 475 END IF 476 477 478 545 FORMAT(//,3X,'Cell Dynamics Parameters (from STDIN)',/ & 479 ,3X,'-------------------------------------') 480 546 FORMAT( 3X,'Simulation cell read from STDIN') 481 547 FORMAT( 3X,'Starting cell generated from CELLDM') 482 548 FORMAT( 3X,'Cell parameters will be re-read from restart file') 483 525 FORMAT( 3X,'Constant VOLUME Molecular dynamics') 484 606 format( 3X,'cell parameters are not allowed to move') 485 526 FORMAT( 3X,'Volume dynamics with steepest descent') 486 527 FORMAT( 3X,'Volume dynamics with newton equations') 487 530 FORMAT( 3X,'Constant PRESSURE Molecular dynamics:',/ & 488 ,3X,'External pressure (GPa) = ',F11.2,/ & 489 ,3X,'Volume mass = ',F11.2) 490 563 FORMAT( 3X,'Zero initial momentum for cell variables') 491 565 FORMAT( 3X,'Volume dynamics: the temperature is not controlled') 492 604 format( 3X,'cell parameters dynamics with nose` temp. control' ) 493 494 600 format( 3X, 'internal stress tensor calculated') 495 602 format( 3X, 'cell parameters dynamics with frich = ',f7.4, & 496 & 3X, 'and greash = ',f7.4 ) 497 498 return 499 end subroutine cell_print_info 500 501 502!---------------------------------------------- 503SUBROUTINE gmeshinfo( ) 504!---------------------------------------------- 505 ! 506 ! Print out the number of g vectors for the different mesh 507 ! 508 USE kinds, ONLY: DP 509 USE mp_global, ONLY: nproc_bgrp, intra_bgrp_comm 510 USE io_global, ONLY: ionode, ionode_id, stdout 511 USE mp, ONLY: mp_max, mp_gather 512 use smallbox_gvec, only: ngb 513 USE gvecw, only: ngw_g, ngw, ngwx 514 USE gvecs, only: ngms_g, ngms, ngsx 515 USE gvect, only: ngm, ngm_g, ngmx 516 USE fft_base, ONLY: dfftp, dffts 517 518 IMPLICIT NONE 519 520 INTEGER :: ip, ng_snd(3), ng_rcv( 3, nproc_bgrp ) 521 INTEGER :: ierr, min_val, max_val, i 522 REAL(DP) :: avg_val 523 524 IF( ngm /= dfftp%ngm ) THEN 525 CALL errore( " gmeshinfo ", " number of G-vectors in module gvect not consistent with FFT descriptor ", 1 ) 526 END IF 527 IF( ngms /= dffts%ngm ) THEN 528 CALL errore( " gmeshinfo ", " number of G-vectors in module gvecs not consistent with FFT descriptor ", 2 ) 529 END IF 530 IF( ngw /= dffts%ngw ) THEN 531 CALL errore( " gmeshinfo ", " number of G-vectors in module gvecw not consistent with FFT descriptor ", 2 ) 532 END IF 533 534 IF(ionode) THEN 535 WRITE( stdout,*) 536 WRITE( stdout,*) ' Reciprocal Space Mesh' 537 WRITE( stdout,*) ' ---------------------' 538 END IF 539 540 ng_snd(1) = ngm_g 541 ng_snd(2) = ngm 542 ng_snd(3) = ngmx 543 CALL mp_gather(ng_snd, ng_rcv, ionode_id, intra_bgrp_comm) 544 ! 545 IF(ionode) THEN 546 min_val = MINVAL( ng_rcv(2,:) ) 547 max_val = MAXVAL( ng_rcv(2,:) ) 548 avg_val = REAL(SUM( ng_rcv(2,:) ))/nproc_bgrp 549 WRITE( stdout,1000) 550 WRITE( stdout,1011) ng_snd(1), min_val, max_val, avg_val 551 END IF 552 ! 553 ng_snd(1) = ngms_g 554 ng_snd(2) = ngms 555 ng_snd(3) = ngsx 556 CALL mp_gather(ng_snd, ng_rcv, ionode_id, intra_bgrp_comm) 557 ! 558 ierr = 0 559 ! 560 IF(ionode) THEN 561 WRITE( stdout,1001) 562 min_val = MINVAL( ng_rcv(2,:) ) 563 max_val = MAXVAL( ng_rcv(2,:) ) 564 avg_val = REAL(SUM( ng_rcv(2,:) ))/nproc_bgrp 565 WRITE( stdout,1011) ng_snd(1), min_val, max_val, avg_val 566 IF( min_val < 1 ) ierr = ip 567 END IF 568 ! 569 CALL mp_max( ierr, intra_bgrp_comm ) 570 ! 571 IF( ierr > 0 ) & 572 CALL errore( " gmeshinfo ", " Wow! some processors have no G-vectors ", ierr ) 573 ! 574 ng_snd(1) = ngw_g 575 ng_snd(2) = ngw 576 ng_snd(3) = ngwx 577 CALL mp_gather(ng_snd, ng_rcv, ionode_id, intra_bgrp_comm) 578 ! 579 IF(ionode) THEN 580 WRITE( stdout,1002) 581 min_val = MINVAL( ng_rcv(2,:) ) 582 max_val = MAXVAL( ng_rcv(2,:) ) 583 avg_val = REAL(SUM( ng_rcv(2,:) ))/nproc_bgrp 584 WRITE( stdout,1011) ng_snd(1), min_val, max_val, avg_val 585 IF( min_val < 1 ) ierr = ip 586 END IF 587 ! 588 CALL mp_max( ierr, intra_bgrp_comm ) 589 ! 590 IF( ierr > 0 ) & 591 CALL errore( " gmeshinfo ", " Wow! some processors have no G-vectors ", ierr ) 592 ! 593 IF(ionode .AND. ngb > 0 ) THEN 594 WRITE( stdout,1050) 595 WRITE( stdout,1060) ngb 596 END IF 597 598 1000 FORMAT(3X,'Large Mesh',/, & 599 ' Global(ngm_g) MinLocal MaxLocal Average') 600 1001 FORMAT(3X,'Smooth Mesh',/, & 601 ' Global(ngms_g) MinLocal MaxLocal Average') 602 1002 FORMAT(3X,'Wave function Mesh',/, & 603 ' Global(ngw_g) MinLocal MaxLocal Average') 604 1011 FORMAT( 3I15, F15.2 ) 605 1050 FORMAT(/,3X,'Small box Mesh') 606 1060 FORMAT( 3X, 'ngb = ', I12, ' not distributed to processors' ) 607 608 RETURN 609 610END SUBROUTINE gmeshinfo 611 612!---------------------------------------------- 613SUBROUTINE constraint_info() 614!---------------------------------------------- 615 USE kinds, ONLY: DP 616 USE constraints_module, ONLY: nconstr, constr_tol, & 617 constr_type, constr, constr_target 618 USE io_global, ONLY: ionode, stdout 619 USE control_flags, ONLY: lconstrain 620 ! 621 IMPLICIT NONE 622 ! 623 INTEGER :: ic 624 ! 625 IF( lconstrain .AND. ionode ) THEN 626 ! 627 WRITE( stdout, 10 ) 628 WRITE( stdout, 20 ) nconstr, constr_tol 629 ! 630 DO ic = 1, nconstr 631 ! 632 IF( constr_type( ic ) == 3 ) THEN 633 ! 634 ! distance 635 ! 636 WRITE( stdout, 30 ) ic 637 WRITE( stdout, 40 ) NINT( constr(1,ic) ), & 638 NINT( constr(2,ic) ), constr_target(ic) 639 ! 640 END IF 641 ! 642 END DO 643 ! 644 END IF 645 ! 64610 FORMAT( 3X, "Using constrained dynamics") 64720 FORMAT( 3X, "number of constrain and tolerance: ", I5, D10.2) 64830 FORMAT( 3X, "constrain ", I5, " type distance ") 64940 FORMAT( 3X, " atoms ", I5, I5, " target dist ", F10.5) 650 ! 651END SUBROUTINE constraint_info 652 653 654SUBROUTINE new_atomind_constraints() 655 ! 656 USE kinds, ONLY: DP 657 USE constraints_module, ONLY: constr 658 ! 659 IMPLICIT NONE 660 ! 661 INTEGER :: ic, ia 662 INTEGER :: iaa 663 REAL(DP) :: aa 664 ! 665 ! Substitute the atom index given in the input file 666 ! with the new atom index, after the sort in the 667 ! atomic coordinates. 668 ! 669 DO ic = 1, SIZE( constr, 2 ) 670 DO ia = 1, SIZE( constr, 1 ) 671 IF( constr( ia, ic ) > 0.0d0 ) THEN 672 iaa = NINT( constr( ia, ic ) ) 673 aa = DBLE( iaa ) 674 constr( ia, ic ) = aa 675 END IF 676 END DO 677 END DO 678 ! 679 RETURN 680 ! 681END SUBROUTINE new_atomind_constraints 682 683 684SUBROUTINE compute_stress_x( stress, detot, h, omega ) 685 USE kinds, ONLY : DP 686 IMPLICIT NONE 687 REAL(DP), INTENT(OUT) :: stress(3,3) 688 REAL(DP), INTENT(IN) :: detot(3,3), h(3,3), omega 689 integer :: i, j 690 do i=1,3 691 do j=1,3 692 stress(i,j)=-1.d0/omega*(detot(i,1)*h(j,1)+ & 693 & detot(i,2)*h(j,2)+detot(i,3)*h(j,3)) 694 enddo 695 enddo 696 return 697END SUBROUTINE compute_stress_x 698!----------------------------------------------------------------------- 699subroutine formf( tfirst, eself ) 700 !----------------------------------------------------------------------- 701 702 !computes (a) the self-energy eself of the ionic pseudocharges; 703 ! (b) the form factors of: (i) pseudopotential (vps), 704 ! (ii) ionic pseudocharge (rhops) 705 ! also calculated the derivative of vps with respect to 706 ! g^2 (dvps) 707 ! 708 USE kinds, ONLY : DP 709 use mp, ONLY : mp_sum 710 use control_flags, ONLY : iprint, tpre, iverbosity 711 use io_global, ONLY : stdout 712 use mp_global, ONLY : intra_bgrp_comm 713 USE fft_base, ONLY : dffts 714 use cell_base, ONLY : omega, tpiba2, tpiba 715 use ions_base, ONLY : rcmax, zv, nsp, na 716 use local_pseudo, ONLY : vps, vps0, rhops, dvps, drhops 717 use atom, ONLY : rgrid 718 use uspp_param, ONLY : upf 719 use pseudo_base, ONLY : compute_rhops, formfn, formfa, compute_eself 720 use pseudopotential, ONLY : tpstab, vps_sp, dvps_sp 721 use splines, ONLY : spline 722 use gvect, ONLY : gstart, gg 723 use constants, ONLY : autoev 724 ! 725 implicit none 726 logical :: tfirst 727 real(DP) :: eself, DeltaV0 728 ! 729 real(DP) :: vpsum, rhopsum 730 integer :: is, ig 731 REAL(DP) :: cost1, xg 732 733 call start_clock( 'formf' ) 734 ! 735 IF( .NOT. ALLOCATED( rgrid ) ) & 736 CALL errore( ' formf ', ' rgrid not allocated ', 1 ) 737 IF( .NOT. ALLOCATED( upf ) ) & 738 CALL errore( ' formf ', ' upf not allocated ', 1 ) 739 ! 740 ! calculation of gaussian selfinteraction 741 ! 742 eself = compute_eself( na, zv, rcmax, nsp ) 743 744 if( tfirst .or. ( iverbosity > 2 ) )then 745 WRITE( stdout, 1200 ) eself 746 endif 747 ! 748 1200 format(/,3x,'formf: eself=',f12.5) 749 ! 750 do is = 1, nsp 751 752 IF( tpstab ) THEN 753 ! 754 ! Use interpolation table, with cubic spline 755 ! 756 cost1 = 1.0d0/omega 757 ! 758 IF( gstart == 2 ) THEN 759 vps (1,is) = vps_sp(is)%y(1) * cost1 760 dvps(1,is) = dvps_sp(is)%y(1) * cost1 761 END IF 762 ! 763 DO ig = gstart, dffts%ngm 764 xg = SQRT( gg(ig) ) * tpiba 765 vps (ig,is) = spline( vps_sp(is), xg ) * cost1 766 dvps(ig,is) = spline( dvps_sp(is), xg ) * cost1 767 END DO 768 ! 769 ELSE 770 771 call formfn( rgrid(is)%r, rgrid(is)%rab, & 772 upf(is)%vloc(1:rgrid(is)%mesh), zv(is), rcmax(is), gg, & 773 omega, tpiba2, rgrid(is)%mesh, dffts%ngm, tpre, & 774 vps(:,is), vps0(is), dvps(:,is) ) 775 776! obsolete BHS form 777! call formfa( vps(:,is), dvps(:,is), rc1(is), rc2(is), wrc1(is), wrc2(is), & 778! rcl(:,is,lloc(is)), al(:,is,lloc(is)), bl(:,is,lloc(is)), & 779! zv(is), rcmax(is), g, omega, tpiba2, dffts%ngm, gstart, tpre ) 780 781 END IF 782 ! 783 ! fourier transform of local pp and gaussian nuclear charge 784 ! 785 call compute_rhops( rhops(:,is), drhops(:,is), zv(is), rcmax(is), gg, & 786 omega, tpiba2, dffts%ngm, tpre ) 787 788 if( tfirst .or. ( iverbosity > 2 ) )then 789 vpsum = SUM( vps( 1:dffts%ngm, is ) ) 790 rhopsum = SUM( rhops( 1:dffts%ngm, is ) ) 791 call mp_sum( vpsum, intra_bgrp_comm ) 792 call mp_sum( rhopsum, intra_bgrp_comm ) 793 WRITE( stdout,1250) (vps(ig,is),rhops(ig,is),ig=1,5) 794 WRITE( stdout,1300) vpsum,rhopsum 795 endif 796 ! 797 end do 798 ! 799 ! ... DeltaV0 is the shift to be applied to eigenvalues 800 ! ... in order to align them to other plane wave codes 801 ! 802 DeltaV0 = 0.0_dp 803 DO is = 1, nsp 804 ! 805 ! ... na(is)/omega is the structure factor at G=0 806 ! 807 DeltaV0 = DeltaV0 + na(is) / omega * vps0(is) 808 END DO 809 ! 810 IF ( tfirst .or. ( iverbosity > 2 ) ) THEN 811 write(6,'(" Delta V(G=0): ",f10.6,"Ry, ",f11.6,"eV")') & 812 deltaV0, deltaV0*autoev 813 END IF 814 ! 815 call stop_clock( 'formf' ) 816 ! 817 1250 format(3x,'formf: vps(g=0)=',f12.7,' rhops(g=0)=',f12.7) 818 1300 format(3x,'formf: sum_g vps(g)=',f12.7,' sum_g rhops(g)=',f12.7) 819 ! 820 return 821end subroutine formf 822! 823!----------------------------------------------------------------------- 824SUBROUTINE newnlinit() 825 !----------------------------------------------------------------------- 826 ! 827 ! ... this routine calculates arrays beta, qq, qgb, rhocb 828 ! ... and derivatives w.r.t. cell parameters dbeta 829 ! ... See also comments in nlinit 830 ! 831 use control_flags, ONLY : tpre 832 use pseudopotential, ONLY : tpstab 833 use cp_interfaces, ONLY : interpolate_beta, interpolate_qradb, compute_qradx, compute_betagx, & 834 exact_beta, check_tables, exact_qradb, build_pstab, build_cctab 835 use betax, only : mmx, refg 836 use kinds, only : dp 837 use io_global, only : ionode, stdout 838 ! 839 IMPLICIT NONE 840 ! 841 LOGICAL :: recompute_table 842 REAL(DP) :: gmax 843 ! 844 ! ... initialization for vanderbilt species 845 ! 846 CALL start_clock( 'newnlinit' ) 847 848 IF( tpstab ) THEN 849 850 recompute_table = tpre .AND. check_tables( gmax ) 851 ! 852 IF ( recompute_table ) THEN 853 854 IF( ionode ) & 855 WRITE( stdout, * ) "newnliinit: recomputing the pseudopotentials tables" 856 !"! 857 858 CALL set_interpolation_table_size( mmx, refg, gmax ) 859 860 CALL compute_qradx( tpre ) 861 862 call compute_betagx( tpre ) 863 864 call build_pstab() 865 ! 866 call build_cctab() 867 868 END IF 869 ! 870 ! initialization that is common to all species 871 ! 872 CALL interpolate_beta( tpre ) 873 ! 874 CALL interpolate_qradb( tpre ) 875 ! 876 ELSE 877 ! 878 ! ... this is mainly for testing 879 ! 880 CALL exact_beta( tpre ) 881 ! 882 CALL exact_qradb( tpre ) 883 ! 884 END IF 885 ! 886 ! ... non-linear core-correction ( rhocb(ig,is) ) 887 ! 888 CALL core_charge_ftr( tpre ) 889 890 CALL stop_clock( 'newnlinit' ) 891 ! 892 RETURN 893 ! 894END SUBROUTINE newnlinit 895! 896!----------------------------------------------------------------------- 897subroutine nlfh_x( stress, bec_bgrp, dbec, lambda, idesc ) 898 !----------------------------------------------------------------------- 899 ! 900 ! contribution to the internal stress tensor due to the constraints 901 ! 902 USE kinds, ONLY : DP 903 use uspp, ONLY : nkb, qq_nt, indv_ijkb0 904 use uspp_param, ONLY : nh, nhm, upf 905 use ions_base, ONLY : nat, ityp 906 use electrons_base, ONLY : nbspx, nbsp, nudx, nspin, nupdwn, iupdwn, ibgrp_g2l 907 use cell_base, ONLY : omega, h 908 use constants, ONLY : pi, fpi, au_gpa 909 use io_global, ONLY : stdout 910 use control_flags, ONLY : iverbosity 911 USE mp, ONLY : mp_sum 912 USE mp_global, ONLY : intra_bgrp_comm, inter_bgrp_comm 913 914! 915 implicit none 916 917 include 'laxlib.fh' 918 919 INTEGER, INTENT(IN) :: idesc(:,:) 920 REAL(DP), INTENT(INOUT) :: stress(3,3) 921 REAL(DP), INTENT(IN) :: bec_bgrp( :, : ), dbec( :, :, :, : ) 922 REAL(DP), INTENT(IN) :: lambda( :, :, : ) 923! 924 INTEGER :: i, j, ii, jj, inl, iv, jv, ia, is, iss, nss, istart 925 INTEGER :: jnl, ir, ic, nr, nc, ibgrp_i, nrcx 926 REAL(DP) :: fpre(3,3), TT, T1, T2 927 ! 928 REAL(DP), ALLOCATABLE :: tmpbec(:,:), tmpdh(:,:), temp(:,:), bec(:,:,:) 929 ! 930 nrcx = MAXVAL( idesc( LAX_DESC_NRCX, : ) ) 931 ! 932 ALLOCATE( bec( nkb, nrcx, nspin ) ) 933 ! 934 DO iss = 1, nspin 935 IF( idesc( LAX_DESC_ACTIVE_NODE, iss ) > 0 ) THEN 936 nss = nupdwn( iss ) 937 istart = iupdwn( iss ) 938 ic = idesc( LAX_DESC_IC, iss ) 939 nc = idesc( LAX_DESC_NC, iss ) 940 DO i=1,nc 941 ibgrp_i = ibgrp_g2l( i+istart-1+ic-1 ) 942 IF( ibgrp_i > 0 ) THEN 943 bec( :, i, iss ) = bec_bgrp( :, ibgrp_i ) 944 ELSE 945 bec( :, i, iss ) = 0.0d0 946 END IF 947 END DO 948 ELSE 949 bec(:,:,iss) = 0.0d0 950 END IF 951 END DO 952 953 CALL mp_sum( bec, inter_bgrp_comm ) 954 ! 955 IF (nspin == 1) THEN 956 IF( ( idesc( LAX_DESC_ACTIVE_NODE, 1 ) > 0 ) ) THEN 957 ALLOCATE ( tmpbec(nhm,nrcx), tmpdh(nrcx,nhm), temp(nrcx,nrcx) ) 958 ENDIF 959 ELSEIF (nspin == 2) THEN 960 IF( ( idesc( LAX_DESC_ACTIVE_NODE, 1 ) > 0 ) .OR. ( idesc( LAX_DESC_ACTIVE_NODE, 2 ) > 0 ) ) THEN 961 ALLOCATE ( tmpbec(nhm,nrcx), tmpdh(nrcx,nhm), temp(nrcx,nrcx) ) 962 END IF 963 ENDIF 964 ! 965 fpre = 0.d0 966 ! 967 do ii=1,3 968 969 do jj=1,3 970 971 do ia=1,nat 972 is = ityp(ia) 973 974 IF( upf(is)%tvanp ) THEN 975 976 do iss = 1, nspin 977 ! 978 istart = iupdwn( iss ) 979 nss = nupdwn( iss ) 980 ! 981 IF( idesc( LAX_DESC_ACTIVE_NODE, iss ) > 0 ) THEN 982 983 ic = idesc( LAX_DESC_IC, iss ) 984 nc = idesc( LAX_DESC_NC, iss ) 985 ir = idesc( LAX_DESC_IR, iss ) 986 nr = idesc( LAX_DESC_NR, iss ) 987 988 tmpbec = 0.d0 989 tmpdh = 0.d0 990! 991 do iv=1,nh(is) 992 do jv=1,nh(is) 993 inl=indv_ijkb0(ia) + jv 994 if(abs(qq_nt(iv,jv,is)).gt.1.e-5) then 995 do i = 1, nc 996 tmpbec(iv,i) = tmpbec(iv,i) + qq_nt(iv,jv,is) * bec( inl, i, iss ) 997 end do 998 endif 999 end do 1000 end do 1001 1002 do iv=1,nh(is) 1003 inl=indv_ijkb0(ia) + iv 1004 do i = 1, nr 1005 tmpdh(i,iv) = dbec( inl, i + (iss-1)*nrcx, ii, jj ) 1006 end do 1007 end do 1008 1009 if(nh(is).gt.0)then 1010 1011 CALL dgemm & 1012 ( 'N', 'N', nr, nc, nh(is), 1.0d0, tmpdh, nrcx, tmpbec, nhm, 0.0d0, temp, nrcx ) 1013 1014 do j = 1, nc 1015 do i = 1, nr 1016 fpre(ii,jj) = fpre(ii,jj) + 2D0 * temp( i, j ) * lambda(i,j,iss) 1017 end do 1018 end do 1019 endif 1020 1021 END IF 1022 ! 1023 end do 1024 ! 1025 END IF 1026 ! 1027 end do 1028 ! 1029 end do 1030 ! 1031 end do 1032 1033 CALL mp_sum( fpre, intra_bgrp_comm ) 1034 1035 do i=1,3 1036 do j=1,3 1037 stress(i,j)=stress(i,j)+ & 1038 (fpre(i,1)*h(j,1)+fpre(i,2)*h(j,2)+fpre(i,3)*h(j,3))/omega 1039 enddo 1040 enddo 1041 1042 IF (allocated(tmpbec)) THEN 1043 DEALLOCATE ( tmpbec, tmpdh, temp ) 1044 END IF 1045 1046 DEALLOCATE( bec ) 1047 1048 1049 IF( iverbosity > 1 ) THEN 1050 WRITE( stdout,*) 1051 WRITE( stdout,*) "constraints contribution to stress" 1052 WRITE( stdout,5555) ((-fpre(i,j),j=1,3),i=1,3) 1053 fpre = MATMUL( fpre, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 1054 WRITE( stdout,5555) ((fpre(i,j),j=1,3),i=1,3) 1055 WRITE( stdout,*) 1056 END IF 1057! 1058 10595555 FORMAT(1x,f12.5,1x,f12.5,1x,f12.5/ & 1060 & 1x,f12.5,1x,f12.5,1x,f12.5/ & 1061 & 1x,f12.5,1x,f12.5,1x,f12.5//) 1062 1063 return 1064end subroutine nlfh_x 1065 1066 1067!----------------------------------------------------------------------- 1068subroutine nlinit 1069 !----------------------------------------------------------------------- 1070 ! 1071 ! this routine allocates and initializes arrays beta, qq, qgb, 1072 ! rhocb, and derivatives w.r.t. cell parameters dbeta 1073 ! 1074 ! beta(ig,l,is) = 4pi/sqrt(omega) y^r(l,q^) 1075 ! int_0^inf dr r^2 j_l(qr) betar(l,is,r) 1076 ! 1077 ! Note that beta(g)_lm,is = (-i)^l*beta(ig,l,is) (?) 1078 ! 1079 ! qq_ij=int_0^r q_ij(r)=omega*qg(g=0) 1080 ! 1081 ! beta and qradb are first calculated on a fixed linear grid in |G| 1082 ! (betax, qradx) then calculated on the box grid by interpolation 1083 ! (this is done in routine newnlinit) 1084 ! 1085 use kinds, ONLY : dp 1086 use control_flags, ONLY : iprint, tpre 1087 use io_global, ONLY : stdout, ionode 1088 use gvecw, ONLY : ngw 1089 use core, ONLY : rhocb, allocate_core 1090 use constants, ONLY : pi, fpi 1091 use ions_base, ONLY : na, nsp 1092 use uspp, ONLY : aainit, beta, qq_nt, dvan, nhtol, nhtolm, indv,& 1093 dbeta 1094 use uspp_param, ONLY : upf, lmaxq, nbetam, lmaxkb, nhm, nh, ish 1095 use atom, ONLY : rgrid 1096 use qgb_mod, ONLY : qgb, dqgb 1097 use smallbox_gvec, ONLY : ngb 1098 use cp_interfaces, ONLY : pseudopotential_indexes, compute_dvan, & 1099 compute_betagx, compute_qradx, build_pstab, build_cctab 1100 USE fft_base, ONLY : dfftp 1101 use pseudopotential, ONLY : tpstab 1102 1103! 1104 implicit none 1105! 1106 integer is, il, l, ir, iv, jv, lm, ind, ltmp, i0 1107 real(dp), allocatable:: fint(:), jl(:), jltmp(:), djl(:), & 1108 & dfint(:) 1109 real(dp) xg, xrg, fac 1110 1111 CALL start_clock( 'nlinit' ) 1112 1113 IF( ionode ) THEN 1114 WRITE( stdout, 100 ) 1115 100 FORMAT( //, & 1116 3X,'Pseudopotentials initialization',/, & 1117 3X,'-------------------------------' ) 1118 END IF 1119 1120 IF( .NOT. ALLOCATED( rgrid ) ) & 1121 CALL errore( ' nlinit ', ' rgrid not allocated ', 1 ) 1122 IF( .NOT. ALLOCATED( upf ) ) & 1123 CALL errore( ' nlinit ', ' upf not allocated ', 1 ) 1124 ! 1125 ! initialize indexes 1126 ! 1127 CALL pseudopotential_indexes( ) 1128 ! 1129 ! initialize array ap 1130 ! 1131 call aainit( lmaxkb + 1 ) 1132 ! 1133 CALL allocate_core( dfftp%nnr, dfftp%ngm, ngb, nsp ) 1134 ! 1135 ! 1136 allocate( beta( ngw, nhm, nsp ) ) 1137 allocate( qgb( ngb, nhm*(nhm+1)/2, nsp ) ) 1138 allocate( qq_nt( nhm, nhm, nsp ) ) 1139 qq_nt (:,:,:) =0.d0 1140 IF (tpre) THEN 1141 allocate( dqgb( ngb, nhm*(nhm+1)/2, nsp, 3, 3 ) ) 1142 allocate( dbeta( ngw, nhm, nsp, 3, 3 ) ) 1143 END IF 1144 ! 1145 ! initialization for vanderbilt species 1146 ! 1147 CALL compute_qradx( tpre ) 1148 ! 1149 ! initialization that is common to all species 1150 ! 1151 WRITE( stdout, fmt="(//,3X,'Common initialization' )" ) 1152 1153 do is = 1, nsp 1154 WRITE( stdout, fmt="(/,3X,'Specie: ',I5)" ) is 1155 ! fac converts ry to hartree 1156 fac=0.5d0 1157 do iv = 1, nh(is) 1158 WRITE( stdout,901) iv, indv(iv,is), nhtol(iv,is) 1159 end do 1160 901 format(2x,i2,' indv= ',i2,' ang. mom= ',i2) 1161 ! 1162 WRITE( stdout,*) 1163 WRITE( stdout,'(20x,a)') ' dion ' 1164 do iv = 1, upf(is)%nbeta 1165 WRITE( stdout,'(8f9.4)') ( fac*upf(is)%dion(iv,jv), jv = 1, upf(is)%nbeta ) 1166 end do 1167 ! 1168 end do 1169 ! 1170 ! calculation of array betagx(ig,iv,is) 1171 ! 1172 call compute_betagx( tpre ) 1173 ! 1174 ! calculate array dvan(iv,jv,is) 1175 ! 1176 call compute_dvan() 1177 ! 1178 IF( tpstab ) THEN 1179 1180 call build_pstab() 1181 ! 1182 call build_cctab() 1183 ! 1184 END IF 1185 ! 1186 ! newnlinit stores qgb and qq, calculates arrays beta rhocb 1187 ! and derivatives wrt cell dbeta 1188 ! 1189 call newnlinit() 1190 1191 CALL stop_clock( 'nlinit' ) 1192 1193 return 1194end subroutine nlinit 1195 1196!------------------------------------------------------------------------- 1197subroutine qvan2b(ngy,iv,jv,is,ylm,qg,qradb) 1198 !-------------------------------------------------------------------------- 1199 ! 1200 ! q(g,l,k) = sum_lm (-i)^l ap(lm,l,k) yr_lm(g^) qrad(g,l,l,k) 1201 ! 1202 USE kinds, ONLY : DP 1203 use control_flags, ONLY : iprint, tpre 1204 use uspp, ONLY : nlx, lpx, lpl, ap, indv, nhtolm 1205 use smallbox_gvec, ONLY : ngb 1206 use uspp_param, ONLY : lmaxq, nbetam 1207 use ions_base, ONLY : nsp 1208! 1209 implicit none 1210 ! 1211 integer, intent(in) :: ngy, iv, jv, is 1212 real(DP), intent(in) :: ylm( ngb, lmaxq*lmaxq ) 1213 real(DP), intent(in) :: qradb( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp ) 1214 complex(DP), intent(out) :: qg( ngb ) 1215! 1216 integer :: ivs, jvs, ijvs, ivl, jvl, i, ii, ij, l, lp, ig 1217 complex(DP) :: sig 1218 ! 1219 ! iv = 1..8 s_1 p_x1 p_z1 p_y1 s_2 p_x2 p_z2 p_y2 1220 ! ivs = 1..4 s_1 s_2 p_1 p_2 1221 ! ivl = 1..4 s p_x p_z p_y 1222 ! 1223 ivs=indv(iv,is) 1224 jvs=indv(jv,is) 1225 if (ivs >= jvs) then 1226 ijvs = ivs*(ivs-1)/2 + jvs 1227 else 1228 ijvs = jvs*(jvs-1)/2 + ivs 1229 end if 1230 ! ijvs is the packed index for (ivs,jvs) 1231 ivl=nhtolm(iv,is) 1232 jvl=nhtolm(jv,is) 1233 if (ivl > nlx .OR. jvl > nlx) & 1234 call errore (' qvan2b ', ' wrong dimensions', MAX(ivl,jvl)) 1235 ! 1236 qg(:) = (0.d0, 0.d0) 1237 ! 1238 ! lpx = max number of allowed y_lm 1239 ! lp = composite lm to indentify them 1240 ! 1241 do i=1,lpx(ivl,jvl) 1242 lp=lpl(ivl,jvl,i) 1243 if (lp > lmaxq*lmaxq) call errore(' qvan2b ',' lp out of bounds ',lp) 1244 ! 1245 ! extraction of angular momentum l from lp: 1246 ! l = int ( sqrt( DBLE(l-1) + epsilon) ) + 1 1247 ! 1248 if (lp == 1) then 1249 l=1 1250 else if ((lp >= 2) .and. (lp <= 4)) then 1251 l=2 1252 else if ((lp >= 5) .and. (lp <= 9)) then 1253 l=3 1254 else if ((lp >= 10).and.(lp <= 16)) then 1255 l=4 1256 else if ((lp >= 17).and.(lp <= 25)) then 1257 l=5 1258 else if ((lp >= 26).and.(lp <= 36)) then 1259 l=6 1260 else if ((lp >= 37).and.(lp <= 49)) then 1261 l=7 1262 else 1263 call errore(' qvan2b ',' not implemented ',lp) 1264 endif 1265 ! 1266 ! sig= (-i)^l 1267 ! 1268 sig=(0.d0,-1.d0)**(l-1) 1269 sig=sig*ap(lp,ivl,jvl) 1270 do ig=1,ngy 1271 qg(ig)=qg(ig)+sig*ylm(ig,lp)*qradb(ig,ijvs,l,is) 1272 end do 1273 end do 1274 1275 return 1276end subroutine qvan2b 1277 1278!------------------------------------------------------------------------- 1279subroutine dqvan2b(ngy,iv,jv,is,ylm,dylm,dqg,dqrad,qradb) 1280 !-------------------------------------------------------------------------- 1281 ! 1282 ! dq(i,j) derivatives wrt to h(i,j) of q(g,l,k) calculated in qvan2b 1283 ! 1284 USE kinds, ONLY : DP 1285 use control_flags, ONLY : iprint, tpre 1286 use uspp, ONLY : nlx, lpx, lpl, ap, indv, nhtolm 1287 use smallbox_gvec, ONLY : ngb 1288 use uspp_param, ONLY : lmaxq, nbetam 1289 use ions_base, ONLY : nsp 1290 1291 implicit none 1292 1293 integer, intent(in) :: ngy, iv, jv, is 1294 REAL(DP), INTENT(IN) :: ylm( ngb, lmaxq*lmaxq ), dylm( ngb, lmaxq*lmaxq, 3, 3 ) 1295 complex(DP), intent(out) :: dqg( ngb, 3, 3 ) 1296 REAL(DP), INTENT(IN) :: dqrad( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp, 3, 3 ) 1297 real(DP), intent(in) :: qradb( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp ) 1298 1299 integer :: ivs, jvs, ijvs, ivl, jvl, i, ii, ij, l, lp, ig 1300 complex(DP) :: sig, z1, z2, zfac 1301 ! 1302 ! 1303 ! iv = 1..8 s_1 p_x1 p_z1 p_y1 s_2 p_x2 p_z2 p_y2 1304 ! ivs = 1..4 s_1 s_2 p_1 p_2 1305 ! ivl = 1..4 s p_x p_z p_y 1306 ! 1307 ivs=indv(iv,is) 1308 jvs=indv(jv,is) 1309 ! 1310 if (ivs >= jvs) then 1311 ijvs = ivs*(ivs-1)/2 + jvs 1312 else 1313 ijvs = jvs*(jvs-1)/2 + ivs 1314 end if 1315 ! 1316 ! ijvs is the packed index for (ivs,jvs) 1317 ! 1318 ivl=nhtolm(iv,is) 1319 jvl=nhtolm(jv,is) 1320 ! 1321 if (ivl > nlx .OR. jvl > nlx) & 1322 call errore (' qvan2 ', ' wrong dimensions (2)', MAX(ivl,jvl)) 1323 ! 1324 dqg(:,:,:) = (0.d0, 0.d0) 1325 1326 ! lpx = max number of allowed y_lm 1327 ! lp = composite lm to indentify them 1328 1329 z1 = 0.0d0 1330 z2 = 0.0d0 1331 do i=1,lpx(ivl,jvl) 1332 lp=lpl(ivl,jvl,i) 1333 if (lp > lmaxq*lmaxq) call errore(' dqvan2b ',' lp out of bounds ',lp) 1334 1335 ! extraction of angular momentum l from lp: 1336 ! l = int ( sqrt( DBLE(l-1) + epsilon) ) + 1 1337 ! 1338 if (lp == 1) then 1339 l=1 1340 else if ((lp >= 2) .and. (lp <= 4)) then 1341 l=2 1342 else if ((lp >= 5) .and. (lp <= 9)) then 1343 l=3 1344 else if ((lp >= 10).and.(lp <= 16)) then 1345 l=4 1346 else if ((lp >= 17).and.(lp <= 25)) then 1347 l=5 1348 else if ((lp >= 26).and.(lp <= 36)) then 1349 l=6 1350 else if ((lp >= 37).and.(lp <= 49)) then 1351 l=7 1352 else 1353 call errore(' qvan2b ',' not implemented ',lp) 1354 endif 1355 ! 1356 ! sig= (-i)^l 1357 ! 1358 sig = (0.0d0,-1.0d0)**(l-1) 1359 ! 1360 sig = sig * ap( lp, ivl, jvl ) 1361 ! 1362 do ij=1,3 1363 do ii=1,3 1364 do ig=1,ngy 1365 zfac = ylm(ig,lp) * dqrad(ig,ijvs,l,is,ii,ij) 1366 zfac = zfac - dylm(ig,lp,ii,ij) * qradb(ig,ijvs,l,is) 1367 dqg(ig,ii,ij) = dqg(ig,ii,ij) + sig * zfac 1368 end do 1369 end do 1370 end do 1371 end do 1372 ! 1373 ! WRITE(6,*) 'DEBUG dqvan2b: ', z1, z2 1374 ! 1375 return 1376end subroutine dqvan2b 1377 1378!----------------------------------------------------------------------- 1379subroutine dylmr2_( nylm, ngy, g, gg, ainv, dylm ) 1380 !----------------------------------------------------------------------- 1381 ! 1382 ! temporary CP interface for PW routine dylmr2 1383 ! dylmr2 calculates d Y_{lm} /d G_ipol 1384 ! dylmr2_ calculates G_ipol \sum_k h^(-1)(jpol,k) (dY_{lm} /dG_k) 1385 ! 1386 USE kinds, ONLY: DP 1387 1388 implicit none 1389 ! 1390 integer, intent(IN) :: nylm, ngy 1391 real(DP), intent(IN) :: g (3, ngy), gg (ngy), ainv(3,3) 1392 real(DP), intent(OUT) :: dylm (ngy, nylm, 3, 3) 1393 ! 1394 integer :: ipol, jpol, lm, ig 1395 real(DP), allocatable :: dylmaux (:,:,:) 1396 ! 1397 allocate ( dylmaux(ngy,nylm,3) ) 1398 ! 1399 dylmaux(:,:,:) = 0.d0 1400 ! 1401 do ipol =1,3 1402 call dylmr2 (nylm, ngy, g, gg, dylmaux(1,1,ipol), ipol) 1403 enddo 1404 ! 1405 do ipol =1,3 1406 do jpol =1,3 1407 do lm=1,nylm 1408 do ig = 1, ngy 1409 dylm (ig,lm,ipol,jpol) = (dylmaux(ig,lm,1) * ainv(jpol,1) + & 1410 dylmaux(ig,lm,2) * ainv(jpol,2) + & 1411 dylmaux(ig,lm,3) * ainv(jpol,3) ) & 1412 * g(ipol,ig) 1413 end do 1414 end do 1415 end do 1416 end do 1417 ! 1418 deallocate ( dylmaux ) 1419 ! 1420 return 1421 ! 1422end subroutine dylmr2_ 1423 1424!----------------------------------------------------------------------- 1425 SUBROUTINE denlcc_x( nnr, nspin, vxcr, sfac, drhocg, dcc ) 1426!----------------------------------------------------------------------- 1427! 1428! derivative of non linear core correction exchange energy wrt cell 1429! parameters h 1430! Output in dcc 1431! 1432 USE kinds, ONLY: DP 1433 USE ions_base, ONLY: nsp 1434 USE gvect, ONLY: gstart, g, gg 1435 USE cell_base, ONLY: omega, ainv, tpiba2 1436 USE mp, ONLY: mp_sum 1437 USE mp_global, ONLY: intra_bgrp_comm 1438 USE uspp_param, ONLY: upf 1439 USE fft_interfaces, ONLY: fwfft 1440 USE fft_base, ONLY: dfftp, dffts 1441 USE fft_helper_subroutines, ONLY: fftx_threed2oned 1442 1443 IMPLICIT NONE 1444 1445 ! input 1446 1447 INTEGER, INTENT(IN) :: nnr, nspin 1448 REAL(DP), INTENT(IN) :: vxcr( :, : ) 1449 COMPLEX(DP), INTENT(IN) :: sfac( :, : ) 1450 REAL(DP), INTENT(IN) :: drhocg( :, : ) 1451 1452 ! output 1453 1454 REAL(DP), INTENT(OUT) :: dcc( :, : ) 1455 1456 ! local 1457 1458 INTEGER :: i, j, ig, is 1459 COMPLEX(DP) :: srhoc 1460 REAL(DP) :: vxcc 1461 ! 1462 COMPLEX(DP), ALLOCATABLE :: vxc( : ), vxg(:) 1463! 1464 dcc = 0.0d0 1465 ! 1466 ALLOCATE( vxc( nnr ) ) 1467 ALLOCATE( vxg( dfftp%ngm ) ) 1468 ! 1469 vxc(:) = vxcr(:,1) 1470 ! 1471 IF( nspin > 1 ) vxc(:) = vxc(:) + vxcr(:,2) 1472 ! 1473 CALL fwfft( 'Rho', vxc, dfftp ) 1474 CALL fftx_threed2oned( dfftp, vxc, vxg ) 1475 ! 1476 DO i=1,3 1477 DO j=1,3 1478 DO ig = gstart, dffts%ngm 1479 srhoc = 0.0d0 1480 DO is = 1, nsp 1481 IF( upf(is)%nlcc ) srhoc = srhoc + sfac( ig, is ) * drhocg( ig, is ) 1482 ENDDO 1483 vxcc = DBLE( CONJG( vxg( ig ) ) * srhoc ) / SQRT( gg(ig) * tpiba2 ) 1484 dcc(i,j) = dcc(i,j) + vxcc * & 1485 & 2.d0 * tpiba2 * g(i,ig) * & 1486 & (g(1,ig)*ainv(j,1) + & 1487 & g(2,ig)*ainv(j,2) + & 1488 & g(3,ig)*ainv(j,3) ) 1489 ENDDO 1490 ENDDO 1491 ENDDO 1492 1493 DEALLOCATE( vxg ) 1494 DEALLOCATE( vxc ) 1495 1496 dcc = dcc * omega 1497 1498 CALL mp_sum( dcc( 1:3, 1:3 ), intra_bgrp_comm ) 1499 1500 RETURN 1501 END SUBROUTINE denlcc_x 1502 1503 1504 1505!----------------------------------------------------------------------- 1506 SUBROUTINE dotcsc_x( eigr, cp, ngw, n ) 1507!----------------------------------------------------------------------- 1508! 1509 USE kinds, ONLY: DP 1510 USE ions_base, ONLY: na, nsp, nat, ityp 1511 USE io_global, ONLY: stdout 1512 USE gvect, ONLY: gstart 1513 USE uspp, ONLY: nkb, qq_nt, indv_ijkb0 1514 USE uspp_param, ONLY: nh, ish, upf 1515 USE mp, ONLY: mp_sum 1516 USE mp_global, ONLY: intra_bgrp_comm, nbgrp, inter_bgrp_comm 1517 USE cp_interfaces, ONLY: nlsm1 1518 USE electrons_base, ONLY: ispin, ispin_bgrp, nbspx_bgrp, ibgrp_g2l, iupdwn, nupdwn, nbspx 1519! 1520 IMPLICIT NONE 1521! 1522 INTEGER, INTENT(IN) :: ngw, n 1523 COMPLEX(DP), INTENT(IN) :: eigr(:,:), cp(:,:) 1524! local variables 1525 REAL(DP) rsum, csc(n) ! automatic array 1526 COMPLEX(DP) temp(ngw) ! automatic array 1527 1528 REAL(DP), ALLOCATABLE:: becp(:,:), cp_tmp(:), becp_tmp(:) 1529 INTEGER i,kmax,nnn,k,ig,is,ia,iv,jv,inl,jnl 1530 INTEGER :: ibgrp_i, ibgrp_k 1531! 1532 ALLOCATE( becp( nkb, nbspx_bgrp ) ) 1533 ALLOCATE( cp_tmp( SIZE( cp, 1 ) ) ) 1534 ALLOCATE( becp_tmp( nkb ) ) 1535! 1536! < beta | phi > is real. only the i lowest: 1537! 1538 1539 CALL nlsm1( nbspx_bgrp, 1, nsp, eigr, cp, becp, 2 ) 1540 1541 nnn = MIN( 12, n ) 1542 1543 DO i = nnn, 1, -1 1544 1545 csc = 0.0d0 1546 1547 ibgrp_i = ibgrp_g2l( i ) 1548 IF( ibgrp_i > 0 ) THEN 1549 cp_tmp = cp( :, ibgrp_i ) 1550 ELSE 1551 cp_tmp = 0.0d0 1552 END IF 1553 1554 CALL mp_sum( cp_tmp, inter_bgrp_comm ) 1555 1556 kmax = i 1557! 1558 DO k=1,kmax 1559 ibgrp_k = ibgrp_g2l( k ) 1560 IF( ibgrp_k > 0 ) THEN 1561 DO ig=1,ngw 1562 temp(ig)=CONJG(cp(ig,ibgrp_k))*cp_tmp(ig) 1563 END DO 1564 csc(k)=2.d0*DBLE(SUM(temp)) 1565 IF (gstart == 2) csc(k)=csc(k)-DBLE(temp(1)) 1566 END IF 1567 END DO 1568 1569 CALL mp_sum( csc( 1:kmax ), intra_bgrp_comm ) 1570 1571 IF( ibgrp_i > 0 ) THEN 1572 becp_tmp = becp( :, ibgrp_i ) 1573 ELSE 1574 becp_tmp = 0.0d0 1575 END IF 1576 1577 CALL mp_sum( becp_tmp, inter_bgrp_comm ) 1578 1579 DO k=1,kmax 1580 rsum=0.d0 1581 ibgrp_k = ibgrp_g2l( k ) 1582 IF( ibgrp_k > 0 ) THEN 1583 DO is=1,nsp 1584 IF( .NOT. upf(is)%tvanp ) CYCLE 1585 DO ia=1,nat 1586 IF( ityp(ia) /= is ) CYCLE 1587 DO iv=1,nh(is) 1588 DO jv=1,nh(is) 1589 inl = indv_ijkb0(ia) + iv 1590 jnl = indv_ijkb0(ia) + jv 1591 rsum = rsum + qq_nt(iv,jv,is)*becp_tmp(inl)*becp(jnl,ibgrp_k) 1592 END DO 1593 END DO 1594 END DO 1595 END DO 1596 END IF 1597 csc(k)=csc(k)+rsum 1598 END DO 1599! 1600 CALL mp_sum( csc( 1:kmax ), inter_bgrp_comm ) 1601 WRITE( stdout,'("dotcsc =",12f18.15)') (csc(k),k=1,i) 1602! 1603 END DO 1604 WRITE( stdout,*) 1605! 1606 DEALLOCATE(becp) 1607 DEALLOCATE(cp_tmp) 1608 DEALLOCATE(becp_tmp) 1609! 1610 RETURN 1611 END SUBROUTINE dotcsc_x 1612 1613 1614! 1615!----------------------------------------------------------------------- 1616 FUNCTION enkin_x( c, f, n ) 1617!----------------------------------------------------------------------- 1618 ! 1619 ! calculation of kinetic energy term 1620 ! 1621 USE kinds, ONLY: DP 1622 USE constants, ONLY: pi, fpi 1623 USE gvecw, ONLY: ngw 1624 USE gvect, ONLY: gstart 1625 USE gvecw, ONLY: g2kin 1626 USE mp, ONLY: mp_sum 1627 USE mp_global, ONLY: intra_bgrp_comm 1628 USE cell_base, ONLY: tpiba2 1629 1630 IMPLICIT NONE 1631 1632 REAL(DP) :: enkin_x 1633 1634 ! input 1635 1636 INTEGER, INTENT(IN) :: n 1637 COMPLEX(DP), INTENT(IN) :: c( :, : ) 1638 REAL(DP), INTENT(IN) :: f( : ) 1639 ! 1640 ! local 1641 1642 INTEGER :: ig, i 1643 REAL(DP) :: sk(n) ! automatic array 1644 ! 1645 DO i=1,n 1646 sk(i)=0.0d0 1647 DO ig=gstart,ngw 1648 sk(i)=sk(i)+DBLE(CONJG(c(ig,i))*c(ig,i))*g2kin(ig) 1649 END DO 1650 END DO 1651 1652 CALL mp_sum( sk(1:n), intra_bgrp_comm ) 1653 1654 enkin_x=0.0d0 1655 DO i=1,n 1656 enkin_x=enkin_x+f(i)*sk(i) 1657 END DO 1658 1659 ! ... reciprocal-space vectors are in units of alat/(2 pi) so a 1660 ! ... multiplicative factor (2 pi/alat)**2 is required 1661 1662 enkin_x = enkin_x * tpiba2 1663! 1664 RETURN 1665 END FUNCTION enkin_x 1666 1667!------------------------------------------------------------------------- 1668 SUBROUTINE nlfl_bgrp_x( bec_bgrp, becdr_bgrp, lambda, idesc, fion ) 1669!----------------------------------------------------------------------- 1670! contribution to fion due to the orthonormality constraint 1671! 1672! 1673 USE kinds, ONLY: DP 1674 USE io_global, ONLY: stdout 1675 USE ions_base, ONLY: na, nsp, nat, ityp 1676 USE uspp, ONLY: nhsa=>nkb, qq_nt, indv_ijkb0 1677 USE uspp_param, ONLY: nhm, nh, upf 1678 USE electrons_base, ONLY: nspin, iupdwn, nupdwn, nbspx_bgrp, ibgrp_g2l, i2gupdwn_bgrp, nbspx, & 1679 iupdwn_bgrp, nupdwn_bgrp 1680 USE constants, ONLY: pi, fpi 1681 USE mp, ONLY: mp_sum 1682 USE mp_global, ONLY: intra_bgrp_comm, inter_bgrp_comm 1683! 1684 IMPLICIT NONE 1685 include 'laxlib.fh' 1686 REAL(DP) :: bec_bgrp(:,:), becdr_bgrp(:,:,:) 1687 REAL(DP), INTENT(IN) :: lambda(:,:,:) 1688 INTEGER, INTENT(IN) :: idesc(:,:) 1689 REAL(DP), INTENT(INOUT) :: fion(:,:) 1690 1691! 1692 INTEGER :: k, is, ia, iv, jv, i, j, inl, isa, iss, nss, istart, ir, ic, nr, nc, ibgrp_i 1693 INTEGER :: n1, n2, m1, m2, nrcx 1694 INTEGER :: nrr(nspin), irr, nrrx 1695 REAL(DP), ALLOCATABLE :: temp(:,:), tmpbec(:,:),tmpdr(:,:), tmplam(:,:,:) 1696 REAL(DP), ALLOCATABLE :: fion_tmp(:,:) 1697 REAL(DP), ALLOCATABLE :: bec(:,:,:) 1698 INTEGER, ALLOCATABLE :: ibgrp_l2g(:,:) 1699 ! 1700 CALL start_clock( 'nlfl' ) 1701 ! 1702 ALLOCATE( fion_tmp( 3, nat ) ) 1703 ! 1704 fion_tmp = 0.0d0 1705 ! 1706 nrcx = MAXVAL( idesc( LAX_DESC_NRCX, : ) ) 1707 ! 1708 1709 ! redistribute bec, becdr according to the ortho subgroup 1710 ! this is required because they are combined with "lambda" matrixes 1711 1712 CALL compute_nrr( nrr ) 1713 nrrx = MAXVAL(nrr) 1714 1715 IF( nrrx > 0 ) THEN 1716 ALLOCATE( tmplam( nrrx, nrcx, nspin ) ) 1717 ALLOCATE( ibgrp_l2g( nrrx, nspin ) ) 1718 END IF 1719 1720 CALL get_local_bec() 1721 CALL get_local_lambda() 1722 1723 ! 1724!$omp parallel default(none), & 1725!$omp shared(nrrx,nhm,nrcx,nsp,na,nspin,nrr,nupdwn,iupdwn,idesc,nh,qq_nt,bec,becdr_bgrp,ibgrp_l2g,tmplam,fion_tmp), & 1726!$omp shared(upf, ityp,nat,indv_ijkb0), & 1727!$omp private(tmpdr,temp,tmpbec,is,k,ia,i,iss,nss,istart,ic,nc,jv,iv,inl,ir,nr) 1728 1729 IF( nrrx > 0 ) THEN 1730 ALLOCATE( tmpdr( nrrx, nhm ) ) 1731 ALLOCATE( temp( nrrx, nrcx ) ) 1732 END IF 1733 ALLOCATE( tmpbec( nhm, nrcx ) ) 1734 1735 DO k=1,3 1736 DO is=1,nsp 1737 IF( .NOT. upf(is)%tvanp ) CYCLE 1738!$omp do 1739 DO ia=1,nat 1740 1741 IF( ityp(ia) /= is ) CYCLE 1742 ! 1743 DO iss = 1, nspin 1744 ! 1745 IF( nrr(iss) == 0 ) CYCLE 1746 ! 1747 nss = nupdwn( iss ) 1748 istart = iupdwn( iss ) 1749 ! 1750 tmpbec = 0.d0 1751 ! 1752 IF( idesc( LAX_DESC_ACTIVE_NODE, iss ) > 0 ) THEN 1753 ! tmpbec distributed by columns 1754 ic = idesc( LAX_DESC_IC, iss ) 1755 nc = idesc( LAX_DESC_NC, iss ) 1756 DO jv=1,nh(is) 1757 inl = indv_ijkb0(ia) + jv 1758 DO iv=1,nh(is) 1759 IF(ABS(qq_nt(iv,jv,is)).GT.1.e-5) THEN 1760 DO i=1,nc 1761 tmpbec(iv,i)=tmpbec(iv,i) + qq_nt(iv,jv,is)*bec(inl,i,iss) 1762 END DO 1763 END IF 1764 END DO 1765 END DO 1766 ! tmpdr distributed by rows 1767 ir = idesc( LAX_DESC_IR, iss ) 1768 nr = idesc( LAX_DESC_NR, iss ) 1769 DO iv=1,nh(is) 1770 inl = indv_ijkb0(ia) + iv 1771 DO i=1,nrr(iss) 1772 tmpdr(i,iv) = becdr_bgrp( inl, ibgrp_l2g(i,iss), k ) 1773 END DO 1774 END DO 1775 END IF 1776 ! 1777 IF( nh(is) > 0 )THEN 1778 IF( idesc( LAX_DESC_ACTIVE_NODE, iss ) > 0 ) THEN 1779 nc = idesc( LAX_DESC_NC, iss ) 1780 CALL dgemm( 'N', 'N', nrr(iss), nc, nh(is), 1.0d0, tmpdr, nrrx, tmpbec, nhm, 0.0d0, temp, nrrx ) 1781 DO j = 1, nc 1782 DO i = 1, nrr(iss) 1783 fion_tmp(k,ia) = fion_tmp(k,ia) + 2D0 * temp( i, j ) * tmplam( i, j, iss ) 1784 END DO 1785 END DO 1786 1787 END IF 1788 ENDIF 1789 END DO 1790 END DO 1791!$omp end do 1792 END DO 1793 END DO 1794 ! 1795 DEALLOCATE( tmpbec ) 1796 ! 1797 IF(ALLOCATED(temp)) DEALLOCATE( temp ) 1798 IF(ALLOCATED(tmpdr)) DEALLOCATE( tmpdr ) 1799 1800!$omp end parallel 1801 1802 DEALLOCATE( bec ) 1803 IF(ALLOCATED(tmplam)) DEALLOCATE( tmplam ) 1804 IF(ALLOCATED(ibgrp_l2g)) DEALLOCATE( ibgrp_l2g ) 1805 ! 1806 CALL mp_sum( fion_tmp, inter_bgrp_comm ) 1807 CALL mp_sum( fion_tmp, intra_bgrp_comm ) 1808 ! 1809 fion = fion + fion_tmp 1810 ! 1811 DEALLOCATE( fion_tmp ) 1812 ! 1813 CALL stop_clock( 'nlfl' ) 1814 ! 1815 RETURN 1816 1817 CONTAINS 1818 1819 SUBROUTINE compute_nrr( nrr ) 1820 INTEGER, INTENT(OUT) :: nrr(:) 1821 nrr = 0 1822 DO iss = 1, nspin 1823 nss = nupdwn( iss ) 1824 istart = iupdwn( iss ) 1825 IF( idesc(LAX_DESC_ACTIVE_NODE, iss ) > 0 ) THEN 1826 ir = idesc( LAX_DESC_IR, iss ) 1827 nr = idesc( LAX_DESC_NR, iss ) 1828 DO i=1,nr 1829 ibgrp_i = ibgrp_g2l( i+istart-1+ir-1 ) 1830 IF( ibgrp_i > 0 ) THEN 1831 nrr(iss) = nrr(iss) + 1 1832 END IF 1833 END DO 1834 END IF 1835 END DO 1836 END SUBROUTINE compute_nrr 1837 1838 SUBROUTINE get_local_bec 1839 ALLOCATE( bec( nhsa, nrcx, nspin ) ) 1840 DO iss = 1, nspin 1841 nss = nupdwn( iss ) 1842 istart = iupdwn( iss ) 1843 IF( idesc(LAX_DESC_ACTIVE_NODE, iss ) > 0 ) THEN 1844 ic = idesc( LAX_DESC_IC, iss ) 1845 nc = idesc( LAX_DESC_NC, iss ) 1846 DO i=1,nc 1847 ibgrp_i = ibgrp_g2l( i+istart-1+ic-1 ) 1848 IF( ibgrp_i > 0 ) THEN 1849 bec( :, i, iss ) = bec_bgrp( :, ibgrp_i ) 1850 ELSE 1851 bec( :, i, iss ) = 0.0d0 1852 END IF 1853 END DO 1854 ELSE 1855 bec(:,:,iss) = 0.0d0 1856 END IF 1857 END DO 1858 CALL mp_sum( bec, inter_bgrp_comm ) 1859 END SUBROUTINE get_local_bec 1860 1861 SUBROUTINE get_local_lambda 1862 DO iss = 1, nspin 1863 nss = nupdwn( iss ) 1864 istart = iupdwn( iss ) 1865 IF( idesc(LAX_DESC_ACTIVE_NODE, iss ) > 0 ) THEN 1866 ir = idesc( LAX_DESC_IR, iss ) 1867 nr = idesc( LAX_DESC_NR, iss ) 1868 irr = 0 1869 DO i=1,nr 1870 ibgrp_i = ibgrp_g2l( i+istart-1+ir-1 ) 1871 IF( ibgrp_i > 0 ) THEN 1872 irr = irr + 1 1873 tmplam(irr,:,iss) = lambda(i,:,iss) 1874 ibgrp_l2g(irr,iss) = ibgrp_i 1875 END IF 1876 END DO 1877 tmplam( irr + 1 : nrrx , :, iss ) = 0.0d0 1878 tmplam( 1 : nrrx , idesc( LAX_DESC_NC, iss ) + 1 : nrcx, iss ) = 0.0d0 1879 END IF 1880 END DO 1881 END SUBROUTINE get_local_lambda 1882 1883 END SUBROUTINE nlfl_bgrp_x 1884! 1885!----------------------------------------------------------------------- 1886 SUBROUTINE pbc(rin,a1,a2,a3,ainv,rout) 1887!----------------------------------------------------------------------- 1888! 1889! brings atoms inside the unit cell 1890! 1891 USE kinds, ONLY: DP 1892 1893 IMPLICIT NONE 1894! input 1895 REAL(DP) rin(3), a1(3),a2(3),a3(3), ainv(3,3) 1896! output 1897 REAL(DP) rout(3) 1898! local 1899 REAL(DP) x,y,z 1900! 1901! bring atomic positions to crystal axis 1902! 1903 x = ainv(1,1)*rin(1)+ainv(1,2)*rin(2)+ainv(1,3)*rin(3) 1904 y = ainv(2,1)*rin(1)+ainv(2,2)*rin(2)+ainv(2,3)*rin(3) 1905 z = ainv(3,1)*rin(1)+ainv(3,2)*rin(2)+ainv(3,3)*rin(3) 1906! 1907! bring x,y,z in the range between -0.5 and 0.5 1908! 1909 x = x - NINT(x) 1910 y = y - NINT(y) 1911 z = z - NINT(z) 1912! 1913! bring atomic positions back in cartesian axis 1914! 1915 rout(1) = x*a1(1)+y*a2(1)+z*a3(1) 1916 rout(2) = x*a1(2)+y*a2(2)+z*a3(2) 1917 rout(3) = x*a1(3)+y*a2(3)+z*a3(3) 1918! 1919 RETURN 1920 END SUBROUTINE pbc 1921 1922! 1923!------------------------------------------------------------------------- 1924 SUBROUTINE prefor_x(eigr,betae) 1925!----------------------------------------------------------------------- 1926! 1927! input : eigr = e^-ig.r_i 1928! output: betae_i,i(g) = (-i)**l beta_i,i(g) e^-ig.r_i 1929! 1930 USE kinds, ONLY : DP 1931 USE ions_base, ONLY : nat, ityp 1932 USE gvecw, ONLY : ngw 1933 USE uspp, ONLY : beta, nhtol, indv_ijkb0 1934 USE uspp_param, ONLY : nh, upf 1935! 1936 IMPLICIT NONE 1937 COMPLEX(DP), INTENT(IN) :: eigr( :, : ) 1938 COMPLEX(DP), INTENT(OUT) :: betae( :, : ) 1939! 1940 INTEGER :: is, iv, ia, inl, ig, isa 1941 COMPLEX(DP) :: ci 1942! 1943 CALL start_clock( 'prefor' ) 1944 DO ia=1,nat 1945 is=ityp(ia) 1946 DO iv=1,nh(is) 1947 ci=(0.0d0,-1.0d0)**nhtol(iv,is) 1948 inl = indv_ijkb0(ia) + iv 1949 DO ig=1,ngw 1950 betae(ig,inl)=ci*beta(ig,iv,is)*eigr(ig,ia) 1951 END DO 1952 END DO 1953 END DO 1954 CALL stop_clock( 'prefor' ) 1955! 1956 RETURN 1957 END SUBROUTINE prefor_x 1958 1959!------------------------------------------------------------------------ 1960 SUBROUTINE collect_bec_x( bec_repl, bec_dist, idesc, nspin ) 1961!------------------------------------------------------------------------ 1962 USE kinds, ONLY : DP 1963 USE mp_global, ONLY : intra_bgrp_comm 1964 USE mp, ONLY : mp_sum 1965 USE io_global, ONLY : stdout 1966 IMPLICIT NONE 1967 include 'laxlib.fh' 1968 REAL(DP), INTENT(OUT) :: bec_repl(:,:) 1969 REAL(DP), INTENT(IN) :: bec_dist(:,:) 1970 INTEGER, INTENT(IN) :: idesc(:,:) 1971 INTEGER, INTENT(IN) :: nspin 1972 INTEGER :: i, ir, n, nrcx, iss 1973 ! 1974 bec_repl = 0.0d0 1975 ! 1976 ! bec is distributed across row processor, the first column is enough 1977 ! 1978 IF( idesc( LAX_DESC_ACTIVE_NODE, 1 ) > 0 .AND. ( idesc( LAX_DESC_MYC, 1 ) == 0 ) ) THEN 1979 ir = idesc( LAX_DESC_IR, 1 ) 1980 DO i = 1, idesc( LAX_DESC_NR, 1 ) 1981 bec_repl( :, i + ir - 1 ) = bec_dist( :, i ) 1982 END DO 1983 IF( nspin == 2 ) THEN 1984 n = idesc( LAX_DESC_N, 1 ) ! number of states with spin==1 ( nupdw(1) ) 1985 nrcx = idesc( LAX_DESC_NRCX, 1 ) ! array elements reserved for each spin ( bec(:,2*nrcx) ) 1986 ir = idesc( LAX_DESC_IR, 2 ) 1987 DO i = 1, idesc( LAX_DESC_NR, 2 ) 1988 bec_repl( :, i + ir - 1 + n ) = bec_dist( :, i + nrcx ) 1989 END DO 1990 END IF 1991 END IF 1992 ! 1993 CALL mp_sum( bec_repl, intra_bgrp_comm ) 1994 ! 1995 RETURN 1996 END SUBROUTINE collect_bec_x 1997 ! 1998!------------------------------------------------------------------------ 1999 SUBROUTINE distribute_bec_x( bec_repl, bec_dist, idesc, nspin ) 2000!------------------------------------------------------------------------ 2001 USE kinds, ONLY : DP 2002 IMPLICIT NONE 2003 include 'laxlib.fh' 2004 REAL(DP), INTENT(IN) :: bec_repl(:,:) 2005 REAL(DP), INTENT(OUT) :: bec_dist(:,:) 2006 INTEGER, INTENT(IN) :: idesc(:,:) 2007 INTEGER, INTENT(IN) :: nspin 2008 INTEGER :: i, ir, n, nrcx 2009 ! 2010 IF( idesc( LAX_DESC_ACTIVE_NODE, 1 ) > 0 ) THEN 2011 ! 2012 bec_dist = 0.0d0 2013 ! 2014 ir = idesc( LAX_DESC_IR, 1 ) 2015 DO i = 1, idesc( LAX_DESC_NR, 1 ) 2016 bec_dist( :, i ) = bec_repl( :, i + ir - 1 ) 2017 END DO 2018 ! 2019 IF( nspin == 2 ) THEN 2020 n = idesc( LAX_DESC_N, 1 ) ! number of states with spin 1 ( nupdw(1) ) 2021 nrcx = idesc( LAX_DESC_NRCX, 1 ) ! array elements reserved for each spin ( bec(:,2*nrcx) ) 2022 ir = idesc( LAX_DESC_IR, 2 ) 2023 DO i = 1, idesc( LAX_DESC_NR, 2 ) 2024 bec_dist( :, i + nrcx ) = bec_repl( :, i + ir - 1 + n ) 2025 END DO 2026 END IF 2027 ! 2028 END IF 2029 RETURN 2030 END SUBROUTINE distribute_bec_x 2031