1! 2! Copyright (C) 2002-2005 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! ... wannier function dynamics and electric field 9! - M.S 10! 11!---------------------------------------------------------------------------- 12MODULE efcalc 13 !---------------------------------------------------------------------------- 14 ! 15 USE kinds, ONLY : DP 16 USE io_global, ONLY : stdout 17 USE wannier_base, ONLY : wf_efield, wf_switch 18 USE wannier_base, ONLY : efx0, efy0, efz0, efx1, efy1, efz1, sw_len 19 ! 20 IMPLICIT NONE 21 ! 22 REAL(DP) :: efx, efy, efz, sw_step 23 REAL(DP), ALLOCATABLE :: xdist(:), ydist(:), zdist(:) 24 ! 25 CONTAINS 26 ! 27 !-------------------------------------------------------------------------- 28 SUBROUTINE clear_nbeg( nbeg ) 29 !-------------------------------------------------------------------------- 30 ! 31 ! ... some more electric field stuff 32 ! - M.S 33 ! 34 INTEGER, INTENT(INOUT) :: nbeg 35 ! 36 ! 37 IF ( wf_efield ) THEN 38 ! 39 IF ( wf_switch ) THEN 40 ! 41 WRITE( stdout, '(/,5X,"!----------------------------------!")' ) 42 WRITE( stdout, '( 5X,"! !")' ) 43 WRITE( stdout, '( 5X,"! ADIABATIC SWITCHING OF THE FIELD !")' ) 44 WRITE( stdout, '( 5X,"! !")' ) 45 WRITE( stdout, '( 5X,"!----------------------------------!",/)' ) 46 ! 47 nbeg=0 48 ! 49 END IF 50 ! 51 END IF 52 ! 53 RETURN 54 ! 55 END SUBROUTINE clear_nbeg 56 ! 57 !-------------------------------------------------------------------------- 58 SUBROUTINE ef_force( fion, ityp, nat, zv ) 59 !-------------------------------------------------------------------------- 60 ! 61 ! ... Electric Feild for ions here 62 ! 63 IMPLICIT NONE 64 ! 65 REAL(DP) :: fion(:,:), zv(:) 66 INTEGER :: ityp(:), nat 67 INTEGER :: ia 68 ! 69 IF ( wf_efield ) THEN 70 ! 71 DO ia = 1, nat 72 ! 73 fion(1,ia) = fion(1,ia) + efx * zv(ityp(ia)) 74 fion(2,ia) = fion(2,ia) + efy * zv(ityp(ia)) 75 fion(3,ia) = fion(3,ia) + efz * zv(ityp(ia)) 76 ! 77 END DO 78 ! 79 END IF 80 ! 81 RETURN 82 ! 83 END SUBROUTINE ef_force 84 ! 85 ! 86 SUBROUTINE deallocate_efcalc() 87 IF( ALLOCATED( xdist ) ) DEALLOCATE( xdist ) 88 IF( ALLOCATED( ydist ) ) DEALLOCATE( ydist ) 89 IF( ALLOCATED( zdist ) ) DEALLOCATE( zdist ) 90 END SUBROUTINE deallocate_efcalc 91 ! 92END MODULE efcalc 93! 94!-------------------------------------------------------------------------- 95MODULE tune 96 !-------------------------------------------------------------------------- 97 ! 98 USE kinds, ONLY : DP 99 ! 100 LOGICAL :: shift 101 INTEGER :: npts, av0, av1, xdir, ydir, zdir, start 102 REAL(DP) :: alpha, b 103 ! 104END MODULE tune 105! 106!-------------------------------------------------------------------------- 107MODULE wannier_module 108 !-------------------------------------------------------------------------- 109 ! 110 ! ... In the presence of an electric field every wannier state feels a 111 ! ... different potantial, which depends on the position of its center. 112 ! ... RHOS is read in as the charge density in subrouting vofrho and 113 ! ... overwritten to be the potential. 114 ! ... -M.S 115 ! 116 USE kinds, ONLY : DP 117 ! 118 IMPLICIT NONE 119 ! 120 SAVE 121 ! 122 LOGICAL :: what1, wann_calc 123 REAL(DP) :: wfx, wfy, wfz, ionx, iony, ionz 124 REAL(DP), ALLOCATABLE :: utwf(:,:) 125 REAL(DP), ALLOCATABLE :: wfc(:,:) 126 REAL(DP), ALLOCATABLE :: rhos1(:,:), rhos2(:,:) 127 COMPLEX(DP), ALLOCATABLE :: rhogdum(:,:) 128 ! 129 CONTAINS 130 ! 131 !------------------------------------------------------------------------ 132 SUBROUTINE allocate_wannier( nbsp, nrxxs, nspin, ng ) 133 !------------------------------------------------------------------------ 134 ! 135 INTEGER, INTENT(in) :: nbsp, nrxxs, nspin, ng 136 ! 137 ALLOCATE( utwf( nbsp, nbsp ) ) 138 ALLOCATE( wfc( 3, nbsp ) ) 139 ALLOCATE( rhos1( nrxxs, nspin) ) 140 ALLOCATE( rhos2( nrxxs, nspin) ) 141 ALLOCATE( rhogdum( ng, nspin ) ) 142 ! 143 RETURN 144 ! 145 END SUBROUTINE allocate_wannier 146 ! 147 !------------------------------------------------------------------------ 148 SUBROUTINE deallocate_wannier() 149 !------------------------------------------------------------------------ 150 ! 151 IF ( ALLOCATED( utwf ) ) DEALLOCATE( utwf ) 152 IF ( ALLOCATED( wfc ) ) DEALLOCATE( wfc ) 153 IF ( ALLOCATED( rhos1 ) ) DEALLOCATE( rhos1 ) 154 IF ( ALLOCATED( rhos2 ) ) DEALLOCATE( rhos2 ) 155 IF ( ALLOCATED( rhogdum ) ) DEALLOCATE( rhogdum ) 156 ! 157 RETURN 158 ! 159 END SUBROUTINE deallocate_wannier 160 ! 161END MODULE wannier_module 162! 163!-------------------------------------------------------------------------- 164MODULE electric_field_module 165 !-------------------------------------------------------------------------- 166 ! 167 ! ... 1 Volt / meter = 1/(5.1412*1.e+11) a.u. 168 ! 169 USE kinds, ONLY : DP 170 ! 171 IMPLICIT NONE 172 ! 173 SAVE 174 ! 175 LOGICAL :: field_tune, ft 176 REAL(DP) :: efe_elec, efe_ion, prefactor, e_tuned(3) 177 REAL(DP) :: tt(3), tt2(3) 178 REAL(DP) :: par, alen, blen, clen, rel1(3), rel2(3) 179 ! 180END MODULE electric_field_module 181! 182!-------------------------------------------------------------------------- 183MODULE wannier_subroutines 184 !-------------------------------------------------------------------------- 185 ! 186 USE kinds, ONLY : DP 187 USE io_global, ONLY : stdout, ionode 188 ! 189 IMPLICIT NONE 190 SAVE 191 ! 192 CONTAINS 193 ! 194 !------------------------------------------------------------------------ 195 SUBROUTINE wannier_startup( ibrav, alat, a1, a2, a3, b1, b2, b3 ) 196 !------------------------------------------------------------------------ 197 ! 198 USE wannier_module, ONLY : utwf 199 USE efcalc, ONLY : wf_efield, efx0, efy0, efz0, & 200 efx1, efy1, efz1, wf_switch, sw_len 201 USE wannier_base, ONLY : calwf, wfsd, wfdt, maxwfdt, nsd, nit, & 202 wf_q, wf_friction, nsteps 203 USE printout_base, ONLY : printout_base_name 204 ! 205 IMPLICIT NONE 206 ! 207 INTEGER :: ibrav 208 REAL(DP) :: a1(3), a2(3), a3(3) 209 REAL(DP) :: b1(3), b2(3), b3(3) 210 REAL(DP) :: alat 211 CHARACTER(LEN=256) :: fname 212 ! 213 INTEGER :: i 214 ! 215 ! ... More Wannier and Field Initialization 216 ! 217 IF (calwf.GT.1) THEN 218 IF (calwf.EQ.3 .AND. ionode ) THEN 219 WRITE( stdout, * ) "------------------------DYNAMICS IN THE WANNIER BASIS--------------------------" 220 WRITE( stdout, * ) " DYNAMICS PARAMETERS " 221 IF (wfsd == 1) THEN 222 WRITE( stdout, 12125) wf_q 223 WRITE( stdout, 12126) wfdt 224 WRITE( stdout, 12124) wf_friction 225 WRITE( stdout, * ) nsteps,"STEPS OF DAMPED MOLECULAR DYNAMICS FOR OPTIMIZATION OF THE SPREAD" 226 ELSE IF (wfsd == 2) THEN 227 WRITE( stdout, 12132) wfdt 228 WRITE( stdout, 12133) maxwfdt 229 WRITE( stdout, * ) nsd,"STEPS OF STEEPEST DESCENT FOR OPTIMIZATION OF THE SPREAD" 230 WRITE( stdout, * ) nit-nsd,"STEPS OF CONJUGATE GRADIENT FOR OPTIMIZATION OF THE SPREAD" 231 ELSE 232 WRITE( stdout, * ) "USING JACOBI ROTATIONS FOR OPTIMIZATION OF THE SPREAD" 233 END IF 234 WRITE( stdout, * ) "AVERAGE WANNIER FUNCTION SPREAD WRITTEN TO FORT.24" 235 fname = printout_base_name( "spr" ) 236 WRITE( stdout, * ) "INDIVIDUAL WANNIER FUNCTION SPREAD WRITTEN TO "//TRIM(fname) 237 fname = printout_base_name( "wfc" ) 238 WRITE( stdout, * ) "WANNIER CENTERS WRITTEN TO "//TRIM(fname) 239 WRITE( stdout, * ) "SOME PERTINENT RUN-TIME INFORMATION WRITTEN TO FORT.27" 240 WRITE( stdout, * ) "-------------------------------------------------------------------------------" 241 WRITE( stdout, * ) 24212124 FORMAT(' DAMPING COEFFICIENT USED FOR WANNIER FUNCTION SPREAD OPTIMIZATION = ',f10.7) 24312125 FORMAT(' FICTITIOUS MASS PARAMETER USED FOR SPREAD OPTIMIZATION = ',f7.1) 24412126 FORMAT(' TIME STEP USED FOR DAMPED DYNAMICS = ',f10.7) 245 ! 24612132 FORMAT(' SMALLEST TIMESTEP IN THE SD / CG DIRECTION FOR SPREAD OPTIMIZATION= ',f10.7) 24712133 FORMAT(' LARGEST TIMESTEP IN THE SD / CG DIRECTION FOR SPREAD OPTIMIZATION = ',f10.7) 248 END IF 249 WRITE( stdout, * ) "wannier_startup IBRAV SELECTED:",ibrav 250 ! 251 CALL recips( a1, a2, a3, b1, b2, b3 ) 252 b1 = b1 * alat 253 b2 = b2 * alat 254 b3 = b3 * alat 255 ! 256 CALL wfunc_init( calwf, b1, b2, b3, ibrav) 257 ! 258 WRITE( stdout, * ) 259 utwf=0.d0 260 DO i=1, SIZE( utwf, 1 ) 261 utwf(i, i)=1.d0 262 END DO 263 END IF 264 IF(wf_efield) THEN 265 266 CALL grid_map 267 268 IF( ionode ) THEN 269 WRITE( stdout, * ) "GRID MAPPING DONE" 270 WRITE( stdout, * ) "DYNAMICS IN THE PRESENCE OF AN EXTERNAL ELECTRIC FIELD" 271 WRITE( stdout, * ) 272 WRITE( stdout, * ) "POLARIZATION CONTRIBUTION OUTPUT TO FORT.28 IN THE FOLLOWING FORMAT" 273 WRITE( stdout, * ) 274 WRITE( stdout, * ) "EFX, EFY, EFZ, ELECTRIC ENTHALPY(ELECTRONIC), ELECTRIC ENTHALPY(IONIC)" 275 WRITE( stdout, * ) 276 WRITE( stdout, '(" E0(x) = ",F10.7)' ) efx0 277 WRITE( stdout, '(" E0(y) = ",F10.7)' ) efy0 278 WRITE( stdout, '(" E0(z) = ",F10.7)' ) efz0 279 WRITE( stdout, '(" E1(x) = ",F10.7)' ) efx1 280 WRITE( stdout, '(" E1(y) = ",F10.7)' ) efy1 281 WRITE( stdout, '(" E1(z) = ",F10.7)' ) efz1 282 ! 283 IF ( wf_switch ) WRITE( stdout, 12127) sw_len 284 ! 285 WRITE( stdout, * ) 286 ! 287 END IF 288 ! 28912127 FORMAT(' FIELD WILL BE TURNED ON ADIBATICALLY OVER ',i5,' STEPS') 290 END IF 291 ! 292 RETURN 293 ! 294 END SUBROUTINE wannier_startup 295 ! 296 !-------------------------------------------------------------------------- 297 SUBROUTINE get_wannier_center( tfirst, cm, bec, eigr, & 298 eigrb, taub, irb, ibrav, b1, b2, b3 ) 299 !-------------------------------------------------------------------------- 300 ! 301 USE efcalc, ONLY: wf_efield 302 USE wannier_base, ONLY: calwf, jwf 303 USE wannier_module, ONLY: what1, wfc, utwf 304 ! 305 IMPLICIT NONE 306 ! 307 LOGICAL, INTENT(in) :: tfirst 308 COMPLEX(DP) :: cm(:,:) 309 REAL(DP) :: bec(:,:) 310 COMPLEX(DP) :: eigrb(:,:), eigr(:,:) 311 INTEGER :: irb(:,:) 312 REAL(DP) :: taub(:,:) 313 INTEGER :: ibrav 314 REAL(DP) :: b1(:), b2(:), b3(:) 315 ! 316 ! ... Get Wannier centers for the first step if wf_efield=true 317 ! 318 IF ( wf_efield ) THEN 319 ! 320 IF ( tfirst ) THEN 321 ! 322 what1 = .TRUE. 323 ! 324 jwf = 1 325 ! 326 CALL wf( calwf,cm, bec, eigr, eigrb, taub, irb, & 327 b1, b2, b3, utwf, what1, wfc, jwf, ibrav ) 328 ! 329 what1 = .FALSE. 330 ! 331 END IF 332 END IF 333 ! 334 RETURN 335 ! 336 END SUBROUTINE get_wannier_center 337 ! 338 !-------------------------------------------------------------------------- 339 SUBROUTINE ef_tune( rhog, tau0 ) 340 !-------------------------------------------------------------------------- 341 ! 342 USE electric_field_module, ONLY: field_tune, e_tuned 343 USE wannier_module, ONLY: rhogdum 344 ! 345 IMPLICIT NONE 346 ! 347 COMPLEX(DP) :: rhog(:,:) 348 REAL(DP) :: tau0(:,:) 349 ! 350 ! ... Tune the Electric field 351 ! 352 IF ( field_tune ) THEN 353 ! 354 rhogdum = rhog 355 ! 356 CALL macroscopic_average( rhogdum, tau0, e_tuned ) 357 ! 358 END IF 359 ! 360 RETURN 361 ! 362 END SUBROUTINE ef_tune 363 ! 364 !-------------------------------------------------------------------------- 365 SUBROUTINE write_charge_and_exit( rhog ) 366 !-------------------------------------------------------------------------- 367 ! 368 USE wannier_base, ONLY : writev 369 ! 370 IMPLICIT NONE 371 ! 372 COMPLEX(DP) :: rhog(:,:) 373 ! 374 ! ... Write chargedensity in g-space 375 ! 376 IF ( writev ) THEN 377 ! 378 CALL write_rho_g( rhog ) 379 ! 380 CALL stop_cp_run() 381 ! 382 END IF 383 ! 384 RETURN 385 ! 386 END SUBROUTINE write_charge_and_exit 387 ! 388 !-------------------------------------------------------------------------- 389 SUBROUTINE wf_options( tfirst, nfi, cm, rhovan, bec, dbec, eigr, eigrb, & 390 taub, irb, ibrav, b1, b2, b3, rhor, drhor, rhog, & 391 drhog ,rhos, enl, ekin ) 392 !-------------------------------------------------------------------------- 393 ! 394 USE efcalc, ONLY : wf_efield 395 USE wannier_base, ONLY : nwf, calwf, jwf, wffort, iplot, iwf 396 USE wannier_module, ONLY : what1, wfc, utwf 397 USE cp_interfaces, ONLY : rhoofr 398 USE dener, ONLY : denl, dekin6 399 ! 400 IMPLICIT NONE 401 ! 402 LOGICAL, INTENT(IN) :: tfirst 403 INTEGER :: nfi 404 COMPLEX(DP) :: cm(:,:) 405 REAL(DP) :: bec(:,:) 406 REAL(DP) :: dbec(:,:,:,:) 407 REAL(DP) :: rhovan(:,:,:) 408 COMPLEX(DP) :: eigrb(:,:), eigr(:,:) 409 INTEGER :: irb(:,:) 410 REAL(DP) :: taub(:,:) 411 INTEGER :: ibrav 412 REAL(DP) :: b1(:), b2(:), b3(:) 413 COMPLEX(DP) :: rhog(:,:) 414 COMPLEX(DP) :: drhog(:,:,:,:) 415 REAL(DP) :: drhor(:,:,:,:), rhor(:,:), rhos(:,:) 416 REAL(DP) :: enl, ekin 417 ! 418 INTEGER :: i, j 419 ! 420 ! 421 ! ... Wannier Function options - M.S 422 ! 423 jwf=1 424 IF (calwf.EQ.1) THEN 425 DO i=1, nwf 426 iwf=iplot(i) 427 j=wffort+i-1 428 CALL rhoofr (nfi,cm, irb, eigrb,bec,dbec,rhovan,rhor,drhor,rhog,drhog,rhos,enl,denl,ekin,dekin6,.false.,j) 429 END DO 430 ! 431 CALL stop_cp_run() 432 ! 433 END IF 434 ! 435 IF ( calwf == 2 ) THEN 436 ! 437 ! ... calculate the overlap matrix 438 ! 439 jwf=1 440 ! 441 CALL wf (calwf,cm,bec,eigr,eigrb,taub,irb,b1,b2,b3,utwf,what1,wfc,jwf,ibrav) 442 ! 443 CALL stop_cp_run( ) 444 ! 445 END IF 446 ! 447 IF (calwf.EQ.5) THEN 448 ! 449 jwf=iplot(1) 450 CALL wf (calwf,cm,bec,eigr,eigrb,taub,irb,b1,b2,b3,utwf,what1,wfc,jwf,ibrav) 451 ! 452 CALL stop_cp_run( ) 453 ! 454 END IF 455 ! 456 ! ... End Wannier Function options - M.S 457 ! 458 RETURN 459 END SUBROUTINE wf_options 460 ! 461 !-------------------------------------------------------------------------- 462 SUBROUTINE ef_potential( nfi, rhos, bec, deeq, betae, c0, cm, emadt2, emaver, verl1, verl2 ) 463 !-------------------------------------------------------------------------- 464 ! 465 USE efcalc, ONLY : wf_efield, efx, efy, efz, & 466 efx0, efy0, efz0, efx1, efy1, efz1, & 467 wf_switch, sw_len, sw_step, xdist, & 468 ydist, zdist 469 USE electric_field_module, ONLY : field_tune, e_tuned, par, rel1, rel2 470 USE wannier_module, ONLY : rhos1, rhos2, wfc 471 USE electrons_base, ONLY : nbsp, nspin, nupdwn, f, ispin 472 USE cell_base, ONLY : ainv, alat, at 473 USE gvect, ONLY : gstart 474 USE control_flags, ONLY : tsde 475 USE wave_base, ONLY : wave_steepest, wave_verlet 476 USE cp_interfaces, ONLY : dforce 477 USE fft_base, ONLY : dffts 478 ! 479 IMPLICIT NONE 480 ! 481 INTEGER :: nfi 482 REAL(DP) :: rhos(:,:) 483 REAL(DP) :: bec(:,:) 484 REAL(DP) :: deeq(:,:,:,:) 485 COMPLEX(DP) :: betae(:,:) 486 COMPLEX(DP) :: c0( :, : ) 487 COMPLEX(DP) :: cm( :, : ) 488 REAL(DP) :: emadt2(:) 489 REAL(DP) :: emaver(:) 490 REAL(DP) :: verl1, verl2 491 REAL(DP) :: a1(3), a2(3), a3(3) 492 COMPLEX(DP), ALLOCATABLE :: c2( : ), c3( : ) 493 INTEGER :: i, ir 494 ! 495 ! ... Potential for electric field 496 ! 497 ALLOCATE( c2( SIZE( c0, 1 ))) 498 ALLOCATE( c3( SIZE( c0, 1 ))) 499 500 a1(:) = at(:,1)/alat ; a2(:) = at(:,2)/alat ; a3(:) = at(:,3)/alat 501 502 IF(wf_efield) THEN 503 IF(field_tune) THEN 504 efx=e_tuned(1) 505 efy=e_tuned(2) 506 efz=e_tuned(3) 507 WRITE( stdout, '(" wf_efield Now ",3(F12.8,1X))' ) efx, efy,efz 508 ! 509 ELSE 510 IF(wf_switch) THEN 511 par=0.d0 512 IF(nfi.LE.sw_len) THEN 513 sw_step=1.0d0/DBLE(sw_len) 514 par=nfi*sw_step 515 IF(efx1.LT.efx0) THEN 516 efx=efx0-(efx0-efx1)*par**5*(70*par**4-315*par**3+540*par**2-420*par+126) 517 ELSE 518 efx=efx0+(efx1-efx0)*par**5*(70*par**4-315*par**3+540*par**2-420*par+126) 519 END IF 520 IF(efy1.LT.efy0) THEN 521 efy=efy0-(efy0-efy1)*par**5*(70*par**4-315*par**3+540*par**2-420*par+126) 522 ELSE 523 efy=efy0+(efy1-efy0)*par**5*(70*par**4-315*par**3+540*par**2-420*par+126) 524 END IF 525 IF(efz1.LT.efz0) THEN 526 efz=efz0-(efz0-efz1)*par**5*(70*par**4-315*par**3+540*par**2-420*par+126) 527 ELSE 528 efz=efz0+(efz1-efz0)*par**5*(70*par**4-315*par**3+540*par**2-420*par+126) 529 END IF 530 END IF 531 ELSE 532 efx=efx1 533 efy=efy1 534 efz=efz1 535 END IF 536 END IF 537 END IF 538 DO i=1,nbsp,2 539 IF(wf_efield) THEN 540 rhos1=0.d0 541 rhos2=0.d0 542 DO ir=1,dffts%nnr 543 rel1(1)=xdist(ir)*a1(1)+ydist(ir)*a2(1)+zdist(ir)*a3(1)-wfc(1,i) 544 rel1(2)=xdist(ir)*a1(2)+ydist(ir)*a2(2)+zdist(ir)*a3(2)-wfc(2,i) 545 rel1(3)=xdist(ir)*a1(3)+ydist(ir)*a2(3)+zdist(ir)*a3(3)-wfc(3,i) 546 ! minimum image convention 547 CALL pbc(rel1,a1,a2,a3,ainv,rel1) 548 IF(nspin.EQ.2) THEN 549 IF(i.LE.nupdwn(1)) THEN 550 rhos1(ir,1)=rhos(ir,1)+efx*rel1(1)+efy*rel1(2)+efz*rel1(3) 551 ELSE 552 rhos1(ir,2)=rhos(ir,2)+efx*rel1(1)+efy*rel1(2)+efz*rel1(3) 553 END IF 554 ELSE 555 rhos1(ir,1)=rhos(ir,1)+efx*rel1(1)+efy*rel1(2)+efz*rel1(3) 556 END IF 557 IF(i.NE.nbsp) THEN 558 rel2(1)=xdist(ir)*a1(1)+ydist(ir)*a2(1)+zdist(ir)*a3(1)-wfc(1,i+1) 559 rel2(2)=xdist(ir)*a1(2)+ydist(ir)*a2(2)+zdist(ir)*a3(2)-wfc(2,i+1) 560 rel2(3)=xdist(ir)*a1(3)+ydist(ir)*a2(3)+zdist(ir)*a3(3)-wfc(3,i+1) 561 ! minimum image convention 562 CALL pbc(rel2,a1,a2,a3,ainv,rel2) 563 IF(nspin.EQ.2) THEN 564 IF(i+1.LE.nupdwn(1)) THEN 565 rhos2(ir,1)=rhos(ir,1)+efx*rel2(1)+efy*rel2(2)+efz*rel2(3) 566 ELSE 567 rhos2(ir,2)=rhos(ir,2)+efx*rel2(1)+efy*rel2(2)+efz*rel2(3) 568 END IF 569 ELSE 570 rhos2(ir,1)=rhos(ir,1)+efx*rel2(1)+efy*rel2(2)+efz*rel2(3) 571 END IF 572 ELSE 573 rhos2(ir,:)=rhos1(ir,:) 574 END IF 575 END DO 576 CALL dforce(i,bec,betae,c0,c2,c3,rhos1,dffts%nnr,ispin,f,nbsp,nspin,rhos2) 577 ELSE 578 CALL dforce(i,bec,betae,c0,c2,c3,rhos,dffts%nnr,ispin,f,nbsp,nspin) 579 END IF 580 IF(tsde) THEN 581 CALL wave_steepest( cm(:, i ), c0(:, i ), emadt2, c2 ) 582 CALL wave_steepest( cm(:, i+1), c0(:, i+1), emadt2, c3 ) 583 ELSE 584 CALL wave_verlet( cm(:, i ), c0(:, i ), verl1, verl2, emaver, c2 ) 585 CALL wave_verlet( cm(:, i+1), c0(:, i+1), verl1, verl2, emaver, c3 ) 586 ENDIF 587 IF (gstart.EQ.2) THEN 588 cm(1, i)=CMPLX(DBLE(cm(1, i)),0.d0,kind=DP) 589 cm(1,i+1)=CMPLX(DBLE(cm(1,i+1)),0.d0,kind=DP) 590 END IF 591 END DO 592 593 DEALLOCATE( c2 ) 594 DEALLOCATE( c3 ) 595 596 RETURN 597 END SUBROUTINE ef_potential 598 ! 599 !-------------------------------------------------------------------- 600 ! ... Electric Field Implementation for Electric Enthalpy 601 ! ... - M.S 602 !-------------------------------------------------------------------- 603 ! 604 !-------------------------------------------------------------------------- 605 SUBROUTINE ef_enthalpy( enthal, tau0 ) 606 !-------------------------------------------------------------------------- 607 ! 608 USE efcalc, ONLY : wf_efield, efx, efy, efz 609 USE electric_field_module, ONLY : efe_elec, efe_ion, tt2, tt 610 USE wannier_module, ONLY : wfx, wfy, wfz, ionx, iony, ionz, wfc 611 USE electrons_base, ONLY : nbsp, f 612 USE cell_base, ONLY : ainv, alat, at 613 USE ions_base, ONLY : nsp, zv, ityp, nat 614 USE io_global, ONLY : ionode 615 ! 616 IMPLICIT NONE 617 ! 618 REAL(DP) :: enthal, tau0(:,:) 619 REAL(DP) :: a1(3), a2(3), a3(3) 620 INTEGER :: i, is, ia 621 ! 622 a1(:) = at(:,1)/alat ; a2(:) = at(:,2)/alat ; a3(:) = at(:,3)/alat 623 IF(wf_efield) THEN 624 ! Electronic Contribution First 625 wfx=0.d0 626 wfy=0.d0 627 wfz=0.d0 628 efe_elec=0.d0 629 DO i=1,nbsp 630 tt2(1)=wfc(1,i) 631 tt2(2)=wfc(2,i) 632 tt2(3)=wfc(3,i) 633 CALL pbc(tt2,a1,a2,a3,ainv,tt2) 634 wfx=wfx+f(i)*tt2(1) 635 wfy=wfy+f(i)*tt2(2) 636 wfz=wfz+f(i)*tt2(3) 637 END DO 638 efe_elec=efe_elec+efx*wfx+efy*wfy+efz*wfz 639 !Then Ionic Contribution 640 ionx=0.d0 641 iony=0.d0 642 ionz=0.d0 643 efe_ion=0.d0 644 DO ia=1,nat 645 is = ityp(ia) 646 tt(1)=tau0(1,ia) 647 tt(2)=tau0(2,ia) 648 tt(3)=tau0(3,ia) 649 CALL pbc(tt,a1,a2,a3,ainv,tt) 650 ionx=ionx+zv(is)*tt(1) 651 iony=iony+zv(is)*tt(2) 652 ionz=ionz+zv(is)*tt(3) 653 END DO 654 efe_ion=efe_ion+efx*ionx+efy*iony+efz*ionz 655 IF( ionode ) THEN 656 WRITE(28,'(f12.9,1x,f12.9,1x,f12.9,1x,f20.15,1x,f20.15)') efx, efy, efz, efe_elec,-efe_ion 657 END IF 658 ELSE 659 efe_elec = 0.0_DP 660 efe_ion = 0.0_DP 661 END IF 662 enthal=enthal+efe_elec-efe_ion 663 664 RETURN 665 END SUBROUTINE ef_enthalpy 666 ! 667 !-------------------------------------------------------------------------- 668 SUBROUTINE wf_closing_options( nfi, c0, cm, bec, eigr, eigrb, taub, & 669 irb, ibrav, b1, b2, b3, taus, tausm, vels, & 670 velsm, acc, lambda, lambdam, idesc, xnhe0, xnhem, & 671 vnhe, xnhp0, xnhpm, vnhp, nhpcl,nhpdim,ekincm,& 672 xnhh0, xnhhm, vnhh, velh, ecut, ecutw, delt, & 673 celldm, fion, tps, mat_z, occ_f, rho ) 674 !-------------------------------------------------------------------------- 675 ! 676 USE efcalc, ONLY : wf_efield 677 USE wannier_base, ONLY : nwf, calwf, jwf, wffort, iplot, iwf 678 USE wannier_module, ONLY : what1, wfc, utwf 679 USE electrons_base, ONLY : nbsp 680 USE gvecw, ONLY : ngw 681 USE control_flags, ONLY : ndw 682 USE cell_base, ONLY : h, hold 683 USE uspp, ONLY : nkbus 684 USE cp_interfaces, ONLY : writefile 685 ! 686 IMPLICIT NONE 687 ! 688 INTEGER :: nfi 689 COMPLEX(DP) :: c0(:,:) 690 COMPLEX(DP) :: cm(:,:) 691 REAL(DP) :: bec(:,:) 692 COMPLEX(DP) :: eigrb(:,:), eigr(:,:) 693 INTEGER :: irb(:,:) 694 REAL(DP) :: taub(:,:) 695 INTEGER :: ibrav 696 REAL(DP) :: b1(:), b2(:), b3(:) 697 REAL(DP) :: taus(:,:), tausm(:,:), vels(:,:), velsm(:,:) 698 REAL(DP) :: acc(:) 699 REAL(DP) :: lambda(:,:,:), lambdam(:,:,:) 700 INTEGER, INTENT(IN) :: idesc(:,:) 701 REAL(DP) :: xnhe0, xnhem, vnhe, xnhp0(:), xnhpm(:), vnhp(:), ekincm 702 INTEGER :: nhpcl, nhpdim 703 REAL(DP) :: velh(:,:) 704 REAL(DP) :: xnhh0(:,:), xnhhm(:,:), vnhh(:,:) 705 REAL(DP) :: ecut, ecutw, delt, celldm(:) 706 REAL(DP) :: fion(:,:), tps 707 REAL(DP) :: mat_z(:,:,:), occ_f(:), rho(:,:) 708 ! 709 CALL start_clock('wf_close_opt') 710 ! 711 ! ... More Wannier Function Options 712 ! 713 IF ( calwf == 4 ) THEN 714 ! 715 jwf = 1 716 ! 717 CALL wf( calwf, c0, bec, eigr, eigrb, taub, irb, & 718 b1, b2, b3, utwf, what1, wfc, jwf, ibrav ) 719 ! 720 IF ( nkbus <= 0 ) THEN 721 ! 722 CALL wf( calwf, cm, bec, eigr, eigrb, taub, irb, & 723 b1, b2, b3, utwf, what1, wfc, jwf, ibrav ) 724 ! 725 ELSE 726 ! 727 cm = c0 728 ! 729 END IF 730 ! 731 CALL writefile( h, hold, nfi, c0, cm, taus, & 732 tausm, vels, velsm,acc, lambda, lambdam, idesc, xnhe0, xnhem, & 733 vnhe, xnhp0, xnhpm, vnhp,nhpcl,nhpdim,ekincm, xnhh0, xnhhm,& 734 vnhh, velh, fion, tps, mat_z, occ_f, rho ) 735 ! 736 CALL stop_clock('wf_close_opt') 737 CALL stop_cp_run( ) 738 ! 739 END IF 740 ! 741 IF ( calwf == 3 ) THEN 742 ! 743 ! ... construct overlap matrix and calculate spreads and do Localization 744 ! 745 jwf = 1 746 ! 747 CALL wf( calwf, c0, bec, eigr, eigrb, taub, irb, & 748 b1, b2, b3, utwf, what1, wfc, jwf, ibrav ) 749 ! 750 CALL stop_clock('wf_close_opt') 751 ! 752 END IF 753 ! 754 RETURN 755 ! 756 END SUBROUTINE wf_closing_options 757 ! 758END MODULE wannier_subroutines 759