[Q-e-developers] Re： About Charge Density

[Jiqiang Li]

Dear Paolo,
       I am afraid the potinit.f90 initializes the potential by calling subroutine "v_of_rho" at very beginning of a calculation. However, I want to calculate the potential by new charge density in every SCF loop, |i> changes with the iterations, but the array psi_m(r)*psi_n(r) is fixed. so I focus on the electron.f90 and modified the subroutine "sum_band".  I am not sure the correct place to add the additional term of charge density. I have some puzzle as mentioned in  the previous email:

      (1)The subroutine symrho is marked as "symmetrize rho(G)", what does the "symmetrize charge density in G space' mean? and what's its consequence? Should I add the term before this symmetrization process or after it?

      (2)I add the code after line 443 in subroutine "sum_band" in QE(v5.1) as follow, ( the band index of |m> and |n> are (nbnd -level_ccmg1) and (nbnd-level_ccmg2), respectively. )

IF (loc_ccmg) THEN
    !
    !
    IF ( nks > 1 ) REWIND( iunigk )
    ik_ccmg = 1
    !
    IF (lsda) current_spin = isk(ik_ccmg)
    npw = ngk(ik_ccmg)
    w1_ccmg = 2.0 / omega
    !
    IF (nks > 1) THEN
        READ( iunigk) igk
        CALL get_buffer(evc, nwordwfc, iunwfc, ik_ccmg)  ! copy evc(1:nwordwfc) from the ik-th record in file with a unit of iunwfc, i.e. wavefucntion
    END IF
    !
    IF (nkb > 0) CALL init_us_2(npw, igk, xk(1, ik_ccmg), vkb)
    !
    IF (noncolin) THEN
    !
    ELSE
         IF ( io_wave_ccmg == 'write' ) THEN
        psic(:) = (0.D0, 0.D0)
        psic(nls(igk(1:npw))) =  evc(1:npw, nbnd-level_ccmg1) + &
                                (0.D0, 1.D0) * evc(1:npw, nbnd -level_ccmg2)
        psic(nlsm(igk(1:npw))) = CONJG(evc(1:npw, nbnd -level_ccmg1) - &
                                (0.D0, 1.D0) * evc(1:npw, nbnd -level_ccmg2))
             CALL invfft('Wave', psic, dffts)
             CALL get_rho_gamma_ccmg(rho_r_ccmg(:,current_spin), dffts % nnr, psic, w1_ccmg)
         END IF
     ENDIF ! noncolinear or colinear
 ENDIF  ! loc_ccmg or not
 !
 IF (loc_ccmg) THEN
     IF (noncolin) THEN
     ! the case of noncolinear has not been supported so far, remaining further work.
     ELSE
          IF ( io_wave_ccmg == 'write' ) THEN
              CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid_ccmg, ierr_ccmg)
              WRITE( file_ccmg,'(I4)') myid_ccmg
              file_ccmg = TRIM(ADJUSTL(file_ccmg))
              file_ccmg = 'wave_ccmg/'//'rho'//TRIM(file_ccmg)//'.dat'
              OPEN(unit = 12306, file = file_ccmg)
              DO ig = 1, size(rho_r_ccmg,1)
                  WRITE(12306,*) rho_r_ccmg(ig,current_spin)
              END DO
              CLOSE(12306)
          END IF
          ! if reading wavefunction.
          IF ( io_wave_ccmg == 'read' ) THEN
              CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid_ccmg, ierr_ccmg)
              WRITE( file_ccmg,'(I4)') myid_ccmg
              file_ccmg = TRIM(ADJUSTL(file_ccmg))
              file_ccmg = 'wave_ccmg/'//'rho'//TRIM(file_ccmg)//'.dat'
              INQUIRE(file = file_ccmg, exist = alive_ccmg )
              IF ( alive_ccmg == .false. ) CALL errore( 'sum_band', 'rho.dat does not exsit!', 1 )
              OPEN(unit = 12306, file = file_ccmg)
              DO ig = 1, size(rho_r_ccmg,1)
                  READ(12306,*) rho_r_ccmg(ig,current_spin)
              END DO
              CLOSE(12306)
               rho%of_r(:, current_spin)=rho%of_r(:, current_spin)+rho_r_ccmg(:,current_spin)
          END IF
      ENDIF ! noncolinear or colinear
  ENDIF  ! loc_ccmg or not
  !

 Is the place correct?

Jiqiang
 
 
From: Paolo Giannozzi
Date: 2015-06-12 04:21
To: Jiqiang Li
Subject: Re: About charge density
On Thu, Jun 11, 2015 at 5:49 PM, Jiqiang Li <jqli14 at fudan.edu.cn> wrote:
      
      (3) Actually I want to realize these process:  First I run a whole ordinary SCF calculation and at the end output the  wavefunction on the disk. Then I run another new SCF calculation, but in this new calculation  a new charge density, rho_new = rho + <m|n> = \sum_i(occupation) { <i|i>} + <m|n>, is treated as input charge density in SCF loop . The term <m|n> is fixed as a constant array  read from disk in every SCF loop. So Where should I add the code to realize rho_new=rho + <m|n> ?
 
In subroutine PW/src/potinit.f90, which either reads from file or generates the starting charge 
density (not the potential, actually: it is computed at the end by calling routine "v_of_rho").
You might also use routine "plugin_scf_potential" to add your preferred term to the charge density.

Paolo


-- 
Paolo Giannozzi, Dept. Chemistry&Physics&Environment,
Univ. Udine, via delle Scienze 208, 33100 Udine, Italy
hone +39-0432-558216, fax +39-0432-558222


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20150614/d20ef1a5/attachment.html>