1! 2! Copyright (C) 2002 CP90 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!==================================================================== 10 SUBROUTINE inner_loop_cold( nfi, tfirst, tlast, eigr, irb, eigrb, & 11 rhor, rhog, rhos, rhoc, ei1, ei2, ei3, & 12 sfac, c0, bec, dbec, firstiter, vpot ) 13!==================================================================== 14 ! 15 ! minimizes the total free energy 16 ! using cold smearing, 17 ! 18 ! 19 20 ! declares modules 21 USE kinds, ONLY: dp 22 USE energies, ONLY: eht, epseu, exc, etot, eself, enl, & 23 ekin, atot, entropy, egrand 24 USE electrons_base, ONLY: f, nspin, nel, iupdwn, nupdwn, nudx, & 25 nelt, nx => nbspx, n => nbsp, ispin , & 26 f_bgrp,nupdwn_bgrp,iupdwn_bgrp 27 28 USE ensemble_dft, ONLY: tens, ninner, ismear, etemp, & 29 ef, z0t, c0diag, becdiag, nrcx, & 30 e0, psihpsi, compute_entropy2, & 31 compute_entropy_der, compute_entropy, & 32 niter_cold_restart, lambda_cold 33 USE smallbox_gvec, ONLY: ngb 34 USE gvecw, ONLY: ngw 35 USE gvect, ONLY: gstart 36 USE ions_base, ONLY: na, nat, nsp 37 USE cell_base, ONLY: omega, alat 38 USE fft_base, ONLY: dfftp, dffts 39 USE local_pseudo, ONLY: vps, rhops 40 USE io_global, ONLY: stdout, ionode, ionode_id 41 USE mp_bands, ONLY: intra_bgrp_comm 42 USE dener 43 USE uspp, ONLY: nhsa=> nkb, betae => vkb, & 44 rhovan => becsum, deeq, nlcc_any 45 USE uspp_param, ONLY: nh 46 USE cg_module, ONLY: ene_ok 47 USE ions_positions, ONLY: tau0 48 USE mp, ONLY: mp_sum,mp_bcast, mp_root_sum 49 50 USE cp_interfaces, ONLY: rhoofr, dforce, protate, vofrho, calbec 51 USE cg_module, ONLY: itercg 52 USE cp_main_variables, ONLY: idesc, drhor, drhog 53 ! 54 IMPLICIT NONE 55 56 include 'laxlib.fh' 57 58!input variables 59 INTEGER :: nfi 60 LOGICAL :: tfirst 61 LOGICAL :: tlast 62 COMPLEX(kind=DP) :: eigr( ngw, nat ) 63 COMPLEX(kind=DP) :: c0( ngw, n ) 64 REAL(kind=DP) :: bec( nhsa, n ) 65 REAL(kind=DP) :: dbec( nhsa, n, 3, 3 ) 66 LOGICAL :: firstiter 67 68 69 INTEGER :: irb( 3, nat ) 70 COMPLEX (kind=DP) :: eigrb( ngb, nat ) 71 REAL(kind=DP) :: rhor( dfftp%nnr, nspin ) 72 REAL(kind=DP) :: vpot( dfftp%nnr, nspin ) 73 COMPLEX(kind=DP) :: rhog( dfftp%ngm, nspin ) 74 REAL(kind=DP) :: rhos( dffts%nnr, nspin ) 75 REAL(kind=DP) :: rhoc( dfftp%nnr ) 76 COMPLEX(kind=DP) :: ei1( dfftp%nr1:dfftp%nr1, nat ) 77 COMPLEX(kind=DP) :: ei2( dfftp%nr2:dfftp%nr2, nat ) 78 COMPLEX(kind=DP) :: ei3( dfftp%nr3:dfftp%nr3, nat ) 79 COMPLEX(kind=DP) :: sfac( dffts%ngm, nsp ) 80 81 82!local variables 83 REAL(kind=DP) :: atot0, atot1, atotl, atotmin 84 REAL(kind=DP), ALLOCATABLE :: fion2(:,:), c0hc0(:,:,:) 85 REAL(kind=DP), ALLOCATABLE :: mtmp(:,:) 86 COMPLEX(kind=DP), ALLOCATABLE :: h0c0(:,:) 87 INTEGER :: niter 88 INTEGER :: i,k, is, nss, istart, ig, iss 89 REAL(kind=DP) :: lambda, lambdap 90 REAL(kind=DP), ALLOCATABLE :: epsi0(:,:) 91 92 INTEGER :: np(2), coor_ip(2), ipr, ipc, nr, nc, ir, ic, ii, jj, root, j 93 INTEGER :: idesc_ip(LAX_DESC_SIZE) 94 INTEGER :: np_rot, me_rot, comm_rot, nrlx, leg_ortho 95 96 CALL start_clock( 'inner_loop') 97 98 allocate(fion2(3,nat)) 99 allocate(c0hc0(nrcx, nrcx, nspin)) 100 allocate(h0c0(ngw,nx)) 101 102 CALL laxlib_getval( leg_ortho = leg_ortho ) 103 104 lambdap=0.3d0!small step for free-energy calculation 105 106 107 ! calculates the initial free energy if necessary 108 IF( .not. ene_ok ) THEN 109 110 ! calculates the overlaps bec between the wavefunctions c0 111 ! and the beta functions 112 CALL calbec( 1, nsp, eigr, c0, bec ) 113 114 ! rotates the wavefunctions c0 and the overlaps bec 115 ! (the occupation matrix f_ij becomes diagonal f_i) 116 nrlx = MAXVAL(idesc(LAX_DESC_NRLX,:)) 117 CALL rotate( nrlx, z0t, c0, bec, c0diag, becdiag ) 118 119 ! calculates the electronic charge density 120 CALL rhoofr( nfi, c0diag, irb, eigrb, becdiag, dbec, rhovan, & 121 rhor, drhor, rhog, drhog, rhos, enl, denl, ekin, dekin6 ) 122 IF(nlcc_any) CALL set_cc( irb, eigrb, rhoc ) 123 124 ! calculates the SCF potential, the total energy 125 ! and the ionic forces 126 vpot = rhor 127 CALL vofrho( nfi, vpot, drhor, rhog, drhog, rhos, rhoc, tfirst, & 128 tlast, ei1, ei2, ei3, irb, eigrb, sfac, & 129 tau0, fion2 ) 130 !entropy value already been calculated 131 132 END IF 133 134 atot0=etot+entropy 135 136!starts inner loop 137 do niter=1,ninner 138!calculates c0hc0, which defines the search line (1-\labda)* psihpsi+\labda*c0hc0 139 140 ! calculateas the energy contribution associated with 141 ! the augmentation charges and the 142 ! corresponding contribution to the ionic force 143 144 CALL newd( vpot, irb, eigrb, rhovan, fion2 ) 145 146 ! operates the Hamiltonian on the wavefunction c0 147 h0c0( :, : )= 0.D0 148 DO i= 1, n, 2 149 CALL dforce( i, bec, betae, c0, h0c0(:,i), h0c0(:,i+1), rhos, dffts%nnr, ispin, f, n, nspin ) 150 END DO 151 152 153 154 ! calculates the Hamiltonian matrix in the basis {c0} 155 c0hc0(:,:,:)=0.d0 156 ! 157 DO is= 1, nspin 158 159 nss= nupdwn( is ) 160 istart= iupdwn( is ) 161 162 np(1) = idesc( LAX_DESC_NPR, is ) 163 np(2) = idesc( LAX_DESC_NPC, is ) 164 165 DO ipc = 1, np(2) 166 DO ipr = 1, np(1) 167 168 coor_ip(1) = ipr - 1 169 coor_ip(2) = ipc - 1 170 CALL laxlib_init_desc( idesc_ip, idesc( LAX_DESC_N, is ), idesc( LAX_DESC_NX, is ), np, coor_ip, & 171 idesc( LAX_DESC_COMM, is ), idesc( LAX_DESC_CNTX, is ), 1 ) 172 173 nr = idesc_ip(LAX_DESC_NR) 174 nc = idesc_ip(LAX_DESC_NC) 175 ir = idesc_ip(LAX_DESC_IR) 176 ic = idesc_ip(LAX_DESC_IC) 177 178 CALL GRID2D_RANK( 'R', idesc_ip(LAX_DESC_NPR), idesc_ip(LAX_DESC_NPC), & 179 idesc_ip(LAX_DESC_MYR), idesc_ip(LAX_DESC_MYC), root ) 180 ! 181 root = root * leg_ortho 182 183 ALLOCATE( mtmp( nr, nc ) ) 184 mtmp = 0.0d0 185 186 CALL dgemm( 'T', 'N', nr, nc, 2*ngw, - 2.0d0, c0( 1, istart + ir - 1 ), 2*ngw, & 187 h0c0( 1, istart + ic - 1 ), 2*ngw, 0.0d0, mtmp, nr ) 188 189 IF (gstart == 2) THEN 190 DO jj = 1, nc 191 DO ii = 1, nr 192 i = ii + ir - 1 193 j = jj + ic - 1 194 mtmp(ii,jj) = mtmp(ii,jj) + DBLE( c0( 1, i + istart - 1 ) ) * DBLE( h0c0( 1, j + istart - 1 ) ) 195 END DO 196 END DO 197 END IF 198 199 CALL mp_root_sum( mtmp, c0hc0(1:nr,1:nc,is), root, intra_bgrp_comm ) 200 201! IF( coor_ip(1) == descla( is )%myr .AND. & 202! coor_ip(2) == descla( is )%myc .AND. descla( is )%active_node > 0 ) THEN 203! c0hc0(1:nr,1:nc,is) = mtmp 204! END IF 205 206 DEALLOCATE( mtmp ) 207 208 END DO 209 END DO 210 END DO 211 212 213 if(mod(itercg,niter_cold_restart) == 0) then 214!calculates free energy at lamda=1. 215 CALL inner_loop_lambda( nfi, tfirst, tlast, eigr, irb, eigrb, & 216 rhor, rhog, rhos, rhoc, ei1, ei2, ei3, & 217 sfac, c0, bec, dbec, firstiter,psihpsi,c0hc0,1.d0,atot1, vpot) 218!calculates free energy at lamda=lambdap 219 220 CALL inner_loop_lambda( nfi, tfirst, tlast, eigr, irb, eigrb, & 221 rhor, rhog, rhos, rhoc, ei1, ei2, ei3, & 222 sfac, c0, bec, dbec, firstiter,psihpsi,c0hc0,lambdap,atotl, vpot) 223!find minimum point lambda 224 225 CALL three_point_min(atot0,atotl,atot1,lambdap,lambda,atotmin) 226 227 else 228 atotl=atot0 229 atot1=atot0 230 lambda=lambda_cold 231 endif 232 233!calculates free energy and rho at lambda 234 235 236 ! calculates the new matrix psihpsi 237 238 DO is= 1, nspin 239 psihpsi(:,:,is) = (1.d0-lambda) * psihpsi(:,:,is) + lambda * c0hc0(:,:,is) 240 END DO 241 242 ! diagonalize and calculates energies 243 244 e0( : )= 0.D0 245 246 CALL inner_loop_diag( c0, bec, psihpsi, z0t, e0 ) 247 248 !calculates fro e0 the new occupation and the entropy 249 250 CALL efermi( nelt, n, etemp, 1, f, ef, e0, entropy, ismear, nspin ) 251 252 do is=1,nspin 253 f_bgrp(iupdwn_bgrp(is):iupdwn_bgrp(is)+nupdwn_bgrp(is)-1)=f(1:nupdwn_bgrp(is)) 254 enddo 255 256 257 !calculates new charge and new energy 258 259 260 ! calculates the electronic charge density 261 CALL rhoofr( nfi, c0diag, irb, eigrb, becdiag, dbec, rhovan, & 262 rhor, drhor, rhog, drhog, rhos, enl, denl, ekin, dekin6 ) 263 IF(nlcc_any) CALL set_cc( irb, eigrb, rhoc ) 264 265 ! calculates the SCF potential, the total energy 266 ! and the ionic forces 267 vpot = rhor 268 CALL vofrho( nfi, vpot, drhor, rhog, drhog, rhos, rhoc, tfirst, & 269 tlast, ei1, ei2, ei3, irb, eigrb, sfac, & 270 tau0, fion2 ) 271 !entropy value already been calculated 272 if(ionode) then 273 write(37,*) niter 274 write(37,*) atot0,atotl,atot1 275 write(37,*) lambda,atotmin,etot+entropy 276 endif 277 atotl=atot0 278 atot0=etot+entropy 279 280 281 if(lambda==0.d0) exit 282 283 284 285 enddo 286 287 288 atot=atot0 289 290!ATTENZIONE 291!the following is of capital importance 292 ene_ok= .TRUE. 293!but it would be better to avoid it 294 295 296 DEALLOCATE(fion2) 297 DEALLOCATE(c0hc0) 298 DEALLOCATE(h0c0) 299 300 CALL stop_clock( 'inner_loop' ) 301 return 302!==================================================================== 303 END SUBROUTINE inner_loop_cold 304!==================================================================== 305 306 307 308 309 310 SUBROUTINE inner_loop_lambda( nfi, tfirst, tlast, eigr, irb, eigrb, & 311 rhor, rhog, rhos, rhoc, ei1, ei2, ei3, & 312 sfac, c0, bec, dbec, firstiter,c0hc0,c1hc1,lambda, & 313 free_energy, vpot ) 314 315!this subroutine for the energy matrix (1-lambda)c0hc0+labda*c1hc1 316!calculates the corresponding free energy 317 318 319 ! declares modules 320 USE kinds, ONLY: dp 321 USE energies, ONLY: eht, epseu, exc, etot, eself, enl, & 322 ekin, atot, entropy, egrand 323 USE electrons_base, ONLY: f, nspin, nel, iupdwn, nupdwn, nudx, & 324 nelt, nx => nbspx, n => nbsp, ispin ,& 325 f_bgrp,nupdwn_bgrp,iupdwn_bgrp 326 327 USE ensemble_dft, ONLY: tens, ninner, ismear, etemp, & 328 c0diag, becdiag, z0t, nrcx, nrlx 329 USE smallbox_gvec, ONLY: ngb 330 USE gvecw, ONLY: ngw 331 USE gvect, ONLY: gstart 332 USE ions_base, ONLY: na, nat, nsp 333 USE cell_base, ONLY: omega, alat 334 USE local_pseudo, ONLY: vps, rhops 335 USE io_global, ONLY: stdout, ionode, ionode_id 336 USE dener 337 USE uspp, ONLY: nhsa=> nkb, betae => vkb, & 338 rhovan => becsum, deeq, nlcc_any 339 USE cg_module, ONLY: ene_ok 340 USE ions_positions, ONLY: tau0 341 USE mp, ONLY: mp_sum,mp_bcast 342 use cp_interfaces, only: rhoofr, dforce, vofrho 343 USE cp_main_variables, ONLY: idesc, drhor, drhog 344 USE fft_base, ONLY: dfftp, dffts 345 346 ! 347 IMPLICIT NONE 348 349!input variables 350 INTEGER :: nfi 351 LOGICAL :: tfirst 352 LOGICAL :: tlast 353 COMPLEX(kind=DP) :: eigr( ngw, nat ) 354 COMPLEX(kind=DP) :: c0( ngw, n ) 355 REAL(kind=DP) :: bec( nhsa, n ) 356 REAL(kind=DP) :: dbec( nhsa, n, 3, 3 ) 357 LOGICAL :: firstiter 358 359 360 INTEGER :: irb( 3, nat ) 361 COMPLEX (kind=DP) :: eigrb( ngb, nat ) 362 REAL(kind=DP) :: rhor( dfftp%nnr, nspin ) 363 REAL(kind=DP) :: vpot( dfftp%nnr, nspin ) 364 COMPLEX(kind=DP) :: rhog( dfftp%ngm, nspin ) 365 REAL(kind=DP) :: rhos( dffts%nnr, nspin ) 366 REAL(kind=DP) :: rhoc( dfftp%nnr ) 367 COMPLEX(kind=DP) :: ei1( dfftp%nr1:dfftp%nr1, nat ) 368 COMPLEX(kind=DP) :: ei2( dfftp%nr2:dfftp%nr2, nat ) 369 COMPLEX(kind=DP) :: ei3( dfftp%nr3:dfftp%nr3, nat ) 370 COMPLEX(kind=DP) :: sfac( dffts%ngm, nsp ) 371 372 REAL(kind=DP), INTENT(in) :: c0hc0(nrcx,nrcx,nspin) 373 REAL(kind=DP), INTENT(in) :: c1hc1(nrcx,nrcx,nspin) 374 REAL(kind=DP), INTENT(in) :: lambda 375 REAL(kind=DP), INTENT(out) :: free_energy 376 377 378!local variables 379 REAL(kind=DP), ALLOCATABLE :: clhcl(:,:,:), fion2(:,:) 380 INTEGER :: i,k, is, nss, istart, ig 381 REAL(kind=DP), ALLOCATABLE :: eaux(:), faux(:), zauxt(:,:,:) 382 REAL(kind=DP) :: entropy_aux, ef_aux 383 384 CALL start_clock( 'inner_lambda') 385 386 allocate(clhcl(nrcx, nrcx, nspin)) 387 allocate(eaux(nx)) 388 allocate(faux(nx)) 389 allocate(zauxt(nrlx,nudx,nspin)) 390 allocate(fion2(3,nat)) 391 392 393!calculates the matrix clhcl 394 395 DO is= 1, nspin 396 clhcl(:,:,is)=(1.d0-lambda)*c0hc0(:,:,is)+lambda*c1hc1(:,:,is) 397 END DO 398 399 CALL inner_loop_diag( c0, bec, clhcl, zauxt, eaux ) 400 401 faux(:)=f(:) 402 403 CALL efermi( nelt, n, etemp, 1, f, ef_aux, eaux, entropy_aux, ismear, nspin ) 404 405 do is=1,nspin 406 f_bgrp(iupdwn_bgrp(is):iupdwn_bgrp(is)+nupdwn_bgrp(is)-1)=f(1:nupdwn_bgrp(is)) 407 enddo 408 409 ! calculates the electronic charge density 410 CALL rhoofr( nfi, c0diag, irb, eigrb, becdiag, dbec, rhovan, & 411 rhor, drhor, rhog, drhog, rhos, enl, denl, ekin, dekin6 ) 412 IF(nlcc_any) CALL set_cc( irb, eigrb, rhoc ) 413 414 ! calculates the SCF potential, the total energy 415 ! and the ionic forces 416 vpot = rhor 417 CALL vofrho( nfi, vpot, drhor, rhog, drhog, rhos, rhoc, tfirst, & 418 tlast, ei1, ei2, ei3, irb, eigrb, sfac, & 419 tau0, fion2 ) 420 !entropy value already been calculated 421 422 423 free_energy=etot+entropy_aux 424 425 f(:)=faux(:) 426 427 deallocate(clhcl) 428 deallocate(faux) 429 deallocate(eaux) 430 deallocate(zauxt) 431 deallocate(fion2) 432 433 CALL stop_clock( 'inner_lambda') 434 435 return 436 437 END SUBROUTINE inner_loop_lambda 438 439 440 SUBROUTINE three_point_min(y0,yl,y1,l,lambda,amin) 441!calculates the estimate for the minimum 442 443 USE kinds, ONLY : DP 444 445 446 implicit none 447 448 REAL(kind=DP), INTENT(in) :: y0,yl,y1, l 449 REAL(kind=DP), INTENT(out) :: lambda, amin 450 451 452 REAL(kind=DP) :: a,b,c, x_min, y_min 453 454! factors for f(x)=ax**2+b*x+c 455 c=y0 456 b=(yl-y0-y1*l**2.d0+y0*l**2.d0)/(l-l**2.d0) 457 a=y1-y0-b 458 459 460 x_min=-b/(2.d0*a) 461 if( x_min <= 1.d0 .and. x_min >= 0.d0) then 462 y_min=a*x_min**2.d0+b*x_min+c 463 if(y_min <= y0 .and. y_min <= y1) then 464 lambda=x_min 465 amin=y_min 466 else 467 if(y0 < y1) then 468 lambda=0.d0 469 amin=y0 470 else 471 lambda=1.d0 472 amin=y1 473 endif 474 endif 475 else 476 if(y0 < y1) then 477 lambda=0.d0 478 amin=y0 479 else 480 lambda=1.d0 481 amin=y1 482 endif 483 endif 484 485 486 return 487 488 END SUBROUTINE three_point_min 489 490 491!==================================================================== 492 SUBROUTINE inner_loop_diag( c0, bec, psihpsi, z0t, e0 ) 493!==================================================================== 494 ! 495 ! minimizes the total free energy 496 ! using cold smearing, 497 ! 498 ! declares modules 499 USE kinds, ONLY: dp 500 USE energies, ONLY: eht, epseu, exc, etot, eself, enl, & 501 ekin, atot, entropy, egrand 502 USE electrons_base, ONLY: f, nspin, nel, iupdwn, nupdwn, nudx, & 503 nelt, nx => nbspx, n => nbsp, ispin 504 505 USE ensemble_dft, ONLY: tens, ninner, ismear, etemp, & 506 ef, c0diag, becdiag, & 507 compute_entropy2, nrlx, nrcx, & 508 compute_entropy_der, compute_entropy, & 509 niter_cold_restart, lambda_cold 510 USE smallbox_gvec, ONLY: ngb 511 USE gvecw, ONLY: ngw 512 USE gvect, & 513 ONLY: gstart 514 USE ions_base, ONLY: na, nat, nsp 515 USE cell_base, ONLY: omega, alat 516 USE local_pseudo, ONLY: vps, rhops 517 USE io_global, ONLY: stdout, ionode, ionode_id 518 USE mp_global, ONLY: intra_bgrp_comm 519 USE dener 520 USE uspp, ONLY: nhsa=> nkb, betae => vkb, & 521 rhovan => becsum, deeq 522 USE uspp_param, ONLY: nh 523 USE cg_module, ONLY: ene_ok 524 USE ions_positions, ONLY: tau0 525 USE mp, ONLY: mp_sum,mp_bcast, mp_root_sum 526 527 USE cp_interfaces, ONLY: rhoofr, dforce, protate 528 USE cg_module, ONLY: itercg 529 USE cp_main_variables, ONLY: idesc 530 ! 531 IMPLICIT NONE 532 533 include 'laxlib.fh' 534 535 COMPLEX(kind=DP) :: c0( ngw, n ) 536 REAL(kind=DP) :: bec( nhsa, n ) 537 REAL(kind=DP) :: psihpsi( nrcx, nrcx, nspin ) 538 REAL(kind=DP) :: z0t( nrlx, nudx, nspin ) 539 REAL(kind=DP) :: e0( nx ) 540 541 INTEGER :: i,k, is, nss, istart, ig 542 REAL(kind=DP), ALLOCATABLE :: epsi0(:,:), dval(:), zaux(:,:), mtmp(:,:) 543 544 INTEGER :: np(2), coor_ip(2), ipr, ipc, nr, nc, ir, ic, ii, jj, root, j 545 INTEGER :: np_rot, me_rot, comm_rot, nrl, kk 546 547 CALL start_clock( 'inner_diag') 548 549 e0( : )= 0.D0 550 551 DO is = 1, nspin 552 553 istart = iupdwn( is ) 554 nss = nupdwn( is ) 555 np_rot = idesc( LAX_DESC_NPR, is ) * idesc( LAX_DESC_NPC, is ) 556 me_rot = idesc( LAX_DESC_MYPE, is ) 557 nrl = idesc( LAX_DESC_NRL, is ) 558 comm_rot = idesc( LAX_DESC_COMM, is ) 559 560 allocate( dval( nx ) ) 561 562 dval = 0.0d0 563 564 IF( idesc( LAX_DESC_ACTIVE_NODE, is ) > 0 ) THEN 565 ! 566 ALLOCATE( epsi0( nrl, nss ), zaux( nrl, nss ) ) 567 568 CALL blk2cyc_redist( nss, epsi0, nrl, nss, psihpsi(:,:,is), SIZE(psihpsi,1), SIZE(psihpsi,2), idesc(:,is) ) 569 570 CALL dspev_drv( 'V', epsi0, nrl, dval, zaux, nrl, nrl, nss, np_rot, me_rot, comm_rot ) 571 ! 572 IF( me_rot /= 0 ) dval = 0.0d0 573 ! 574 ELSE 575 576 ALLOCATE( epsi0( 1, 1 ), zaux( 1, 1 ) ) 577 578 END IF 579 580 CALL mp_sum( dval, intra_bgrp_comm ) 581 582 DO i = 1, nss 583 e0( i+istart-1 )= dval( i ) 584 END DO 585 z0t(:,:,is) = 0.0d0 586 587 IF( idesc( LAX_DESC_ACTIVE_NODE, is ) > 0 ) THEN 588 !NB zaux is transposed 589 !ALLOCATE( mtmp( nudx, nudx ) ) 590 z0t( 1:nrl , 1:nss, is ) = zaux( 1:nrl, 1:nss ) 591 END IF 592 593 DEALLOCATE( epsi0 , zaux, dval ) 594 595 END DO 596 597 ! rotates the wavefunctions c0 and the overlaps bec 598 ! (the occupation matrix f_ij becomes diagonal f_i) 599 600 CALL rotate ( nrlx, z0t, c0, bec, c0diag, becdiag ) 601 602 CALL stop_clock( 'inner_diag') 603 604 RETURN 605END SUBROUTINE 606