LCOV - code coverage report
Current view: top level - mix - potdis.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 153 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 0.0 %

          Line data    Source code
       1             : MODULE m_potdis
       2             : CONTAINS
       3           0 :   SUBROUTINE potdis(stars,vacuum,atoms,sphhar, input,cell,sym)
       4             :     !
       5             :     !     *****************************************************
       6             :     !     calculates the root-mean-square distance between input
       7             :     !     and output potentials (averaged over unit cell volume)
       8             :     !                                 based on code by   c.l.fu
       9             :     !     *****************************************************
      10             :     !
      11             :     USE m_intgr, ONLY : intgr3, intgz0
      12             :     USE m_constants, ONLY : fpi_const
      13             :     USE m_loddop
      14             :     USE m_cfft
      15             :     USE m_types
      16             :     IMPLICIT NONE
      17             :     TYPE(t_input),INTENT(IN)   :: input
      18             :     TYPE(t_vacuum),INTENT(IN)  :: vacuum
      19             :     TYPE(t_sym),INTENT(IN)     :: sym
      20             :     TYPE(t_stars),INTENT(IN)   :: stars
      21             :     TYPE(t_cell),INTENT(IN)    :: cell
      22             :     TYPE(t_sphhar),INTENT(IN)  :: sphhar
      23             :     TYPE(t_atoms),INTENT(IN)   :: atoms
      24             :     !     .. Local Scalars ..
      25             :     REAL fact,facv,rhs,sumis,sumz
      26             :     COMPLEX phase
      27             :     INTEGER i,i1,i2,i3,id2,id3,io,ip,iter,ivac,j,k1,k2,k3,lh,n,&
      28             :                  nk12,npz,nt,num,na
      29             :     LOGICAL tail
      30             :     !     ..
      31             :     !     .. Local Arrays ..
      32           0 :     COMPLEX rhpw(stars%ng3,2,2),rhv1(vacuum%nmzxyd,stars%ng2-1,2,2,2)
      33           0 :     REAL rhsp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,2,2),rhv0(vacuum%nmzd,2,2,2)
      34           0 :     REAL dis(2),disz(vacuum%nmzd),rh(atoms%jmtd)
      35           0 :     REAL af3((2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1),input%jspins),&
      36           0 :          bf3((2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1),input%jspins), &
      37           0 :          bf2((2*stars%mx1+1)*(2*stars%mx2+1),input%jspins),&
      38           0 :          af2((2*stars%mx1+1)*(2*stars%mx2+1),input%jspins)
      39           0 :     REAL pdis(0:4,0:atoms%ntype,2)
      40             :     !     ..
      41             :     !
      42             : 
      43           0 :     dis = 0.
      44           0 :     pdis=0.0
      45             : 
      46           0 :     tail = .TRUE.
      47           0 :     npz = vacuum%nmz + 1
      48           0 :     nt = (2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1)
      49           0 :     nk12 = (2*stars%mx1+1)*(2*stars%mx2+1)
      50           0 :     fact = cell%omtil/REAL(nt)
      51           0 :     facv = 1.0
      52             :     !     ---> reload potentials
      53           0 :     OPEN (9,file='nrp',form='unformatted',status='unknown')
      54           0 :     REWIND 9
      55           0 :     DO  io = 1,2
      56             :        CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
      57           0 :             9, iter,rhsp(:,0:,:,:,io),rhpw(:,:,io), rhv0(:,:,:,io),rhv1(:,:,:,:,io))
      58             :     ENDDO
      59           0 :     CLOSE (9)
      60           0 :     IF (input%jspins.EQ.1) THEN
      61             :        !       ---> total potential difference
      62           0 :        CALL priv_cdndif(1,1,1,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
      63             :     ELSE
      64             :        !       ---> total input potential
      65           0 :        CALL priv_cdndif(1,1,2,1,1,1,1.0,1.0,input%film, rhsp,rhpw,rhv0,rhv1)
      66             :        !       ---> total output potential
      67           0 :        CALL priv_cdndif(1,1,2,2,2,2,1.0,1.0,input%film, rhsp,rhpw,rhv0,rhv1)
      68             :        !       ---> input 'spin' potential
      69           0 :        CALL priv_cdndif(2,1,2,1,1,1,1.0,-2.0,input%film, rhsp,rhpw,rhv0,rhv1)
      70             :        !       ---> output 'spin' potential
      71           0 :        CALL priv_cdndif(2,1,2,2,2,2,1.0,-2.0,input%film, rhsp,rhpw,rhv0,rhv1)
      72             :        !       ---> potential difference
      73           0 :        CALL priv_cdndif(1,1,1,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
      74             :        !       ---> spin potential difference
      75           0 :        CALL priv_cdndif(2,2,2,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
      76             :     END IF
      77           0 :     DO  num = 1,input%jspins
      78             :        !     ----> m.t. part
      79           0 :        na = 1
      80           0 :        DO n = 1,atoms%ntype
      81           0 :           rh(:atoms%jri(n))=rhsp(:atoms%jri(n),0,n,num,1)*rhsp(:atoms%jri(n),0,n,num,1) *fpi_const
      82           0 :           CALL intgr3(rh,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),pdis(0,n,num))
      83           0 :           pdis(4,n,num)=pdis(0,n,num)
      84           0 :           DO lh = 1,sphhar%nlh(atoms%ntypsy(na))
      85             :              rh(:atoms%jri(n))=rh(:atoms%jri(n))+rhsp(:atoms%jri(n),lh,n,num,1)*&
      86           0 :                   rhsp(:atoms%jri(n),lh,n,num,1)* atoms%rmsh(:atoms%jri(n),n)*atoms%rmsh(:atoms%jri(n),n)
      87           0 :              CALL intgr3(rh,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),rhs)
      88           0 :              IF (lh<4) pdis(lh,n,num)=rhs
      89           0 :              pdis(4,n,num)=pdis(4,n,num)+rhs
      90             :           ENDDO
      91           0 :           pdis(:,n,num) = pdis(:,n,num)*atoms%neq(n)
      92           0 :           dis(num)=dis(num)+SUM(pdis(:,n,num))
      93           0 :           pdis(:,n,num) = SQRT(pdis(:,n,num)/atoms%volmts(n))*1000.
      94           0 :           na = na + atoms%neq(n)
      95             :        ENDDO
      96             :        !     ----> interstitial part
      97             :        !         ---> create density in the real space
      98           0 :        i = 0
      99           0 :        DO  i3 = 0,2*stars%mx3
     100           0 :           k3 = i3
     101           0 :           IF (k3.GT.stars%mx3) k3 = k3 - 2*stars%mx3-1
     102           0 :           DO  i2 = 0,2*stars%mx2
     103           0 :              k2 = i2
     104           0 :              IF (k2.GT.stars%mx2) k2 = k2 - 2*stars%mx2-1
     105           0 :              DO  i1 = 0,2*stars%mx1
     106           0 :                 k1 = i1
     107           0 :                 IF (k1.GT.stars%mx1) k1 = k1 - 2*stars%mx3-1
     108           0 :                 i = i + 1
     109           0 :                 id3 = stars%ig(k1,k2,k3)
     110           0 :                 phase = stars%rgphs(k1,k2,k3)
     111           0 :                 IF (id3.EQ.0) THEN
     112           0 :                    af3(i,1) = 0.
     113           0 :                    bf3(i,1) = 0.
     114             :                 ELSE
     115             :                    af3(i,1) =  REAL(rhpw(id3,num,1))*REAL(phase) &
     116           0 :                            - AIMAG(rhpw(id3,num,1))*AIMAG(phase)
     117             :                    bf3(i,1) = AIMAG(rhpw(id3,num,1))*REAL(phase) &
     118           0 :                            +  REAL(rhpw(id3,num,1))*AIMAG(phase)
     119             :                 END IF
     120             :              ENDDO
     121             :           ENDDO
     122             :        ENDDO
     123             : 
     124           0 :        CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx1*2+1,stars%mx1*2+1,1)
     125           0 :        CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx2*2+1,nk12,1)
     126           0 :        CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx3*2+1,nt,1)
     127             :        !         ---> form the dot product
     128           0 :        i = 0
     129           0 :        DO  i3 = 0,stars%mx3*2
     130           0 :           DO  i2 = 0,stars%mx2*2
     131           0 :              DO  i1 = 0,stars%mx3*2
     132           0 :                 i = i + 1
     133           0 :                 af3(i,1) = af3(i,1)*af3(i,1)
     134           0 :                 bf3(i,1) = 0.
     135             :              ENDDO
     136             :           ENDDO
     137             :        ENDDO
     138             :        !         ---> back fft
     139           0 :        CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx1*2+1,stars%mx1*2+1,-1)
     140           0 :        CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx2*2+1,nk12,-1)
     141           0 :        CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx3*2+1,nt,-1)
     142           0 :        sumis = 0.
     143           0 :        i = 0
     144           0 :        DO  i3 = 0,stars%mx3*2
     145           0 :           k3 = i3
     146           0 :           IF (k3.GT.stars%mx3) k3 = k3 - stars%mx3*2-1
     147           0 :           DO  i2 = 0,stars%mx2*2
     148           0 :              k2 = i2
     149           0 :              IF (k2.GT.stars%mx2) k2 = k2 - stars%mx2*2-1
     150           0 :              DO  i1 = 0,stars%mx1*2
     151           0 :                 k1 = i1
     152           0 :                 IF (k1.GT.stars%mx1) k1 = k1 - stars%mx1*2-1
     153           0 :                 i = i + 1
     154           0 :                 phase = stars%rgphs(k1,k2,k3)
     155           0 :                 id3 = stars%ig(k1,k2,k3)
     156           0 :                 IF (id3.NE.0) THEN
     157           0 :                    sumis = sumis + REAL(CMPLX(af3(i,1),bf3(i,1))* CONJG(stars%ustep(id3)))*phase*fact
     158             :                 END IF
     159             :              ENDDO
     160             :           ENDDO
     161             :        ENDDO
     162           0 :        pdis(0,0,num)=SQRT(sumis/cell%volint)*1000.
     163           0 :        dis(num) = dis(num) + sumis
     164           0 :        IF (input%film) THEN
     165             :           !     ----> vacuum part
     166           0 :           DO  ivac = 1,vacuum%nvac
     167           0 :              DO  ip = 1,vacuum%nmzxy
     168             :                 !         ---> create density in the real space
     169           0 :                 i = 0
     170           0 :                 DO  i2 = 0,stars%mx2*2
     171           0 :                    k2 = i2
     172           0 :                    IF (k2.GT.stars%mx2) k2 = k2 - stars%mx2*2-1
     173           0 :                    DO  i1 = 0,stars%mx1*2
     174           0 :                       k1 = i1
     175           0 :                       IF (k1.GT.stars%mx1) k1 = k1 - stars%mx1*2-1
     176           0 :                       i = i + 1
     177           0 :                       id3 = stars%ig(k1,k2,0)
     178           0 :                       IF (id3.EQ.0) THEN
     179           0 :                          af2(i,1) = 0.
     180           0 :                          bf2(i,1) = 0.
     181             :                       ELSE
     182           0 :                          id2 = stars%ig2(id3)
     183           0 :                          phase = stars%rgphs(k1,k2,0)
     184           0 :                          IF (id2.EQ.1) THEN
     185           0 :                             af2(i,1) = rhv0(ip,ivac,num,1)
     186           0 :                             bf2(i,1) = 0.
     187             :                          ELSE
     188           0 :                             af2(i,1) = REAL(rhv1(ip,id2-1,ivac,num,1))
     189           0 :                             bf2(i,1) =AIMAG(rhv1(ip,id2-1,ivac,num,1))
     190             :                          END IF
     191             :                       END IF
     192             :                    ENDDO
     193             :                 ENDDO
     194           0 :                 CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx1*2+1,stars%mx1*2+1,1)
     195           0 :                 CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx2*2+1,nk12,1)
     196             :                 !         ---> form dot product
     197           0 :                 i = 0
     198           0 :                 DO  i2 = 0,stars%mx2*2
     199           0 :                    DO  i1 = 0,stars%mx1*2
     200           0 :                       i = i + 1
     201           0 :                       af2(i,1) = af2(i,1)*af2(i,1)
     202           0 :                       bf2(i,1) = 0.
     203             :                    ENDDO
     204             :                 ENDDO
     205             :                 !         ---> back fft
     206           0 :                 CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx1*2+1,stars%mx1*2+1,-1)
     207           0 :                 CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx2*2+1,nk12,-1)
     208           0 :                 disz(npz-ip) = cell%area*af2(1,1)/REAL(nk12)
     209             :              ENDDO
     210             :              !         ---> beyond warping region
     211           0 :              DO  ip = vacuum%nmzxy + 1,vacuum%nmz
     212           0 :                 disz(npz-ip) = cell%area*rhv0(ip,ivac,num,1)
     213             :              ENDDO
     214           0 :              CALL intgz0(disz,vacuum%delz,vacuum%nmz,sumz,tail)
     215           0 :              IF (sym%zrfs .OR. sym%invs) facv = 2.0
     216           0 :              dis(num) = dis(num) + facv*sumz
     217             :           ENDDO
     218             :        END IF
     219             :     ENDDO
     220             : 
     221           0 :     dis(:) = SQRT(dis(:)/cell%vol)*1000.
     222             : 
     223           0 :     WRITE (6,FMT=8000) iter,dis(1)
     224             : 8000 FORMAT (/,'----> distance of  the potential for it=',i3,':',f11.6, ' mhtr/bohr**3')
     225           0 :     WRITE(6,*) "Details of potential differences for each atom type"
     226           0 :     WRITE(6,*) "Atom: total difference: difference of first sphhar"
     227           0 :     DO n=1,atoms%ntype
     228           0 :        WRITE(6,"(i5,' : ',f10.6,' : ',4f10.6)") n,pdis(4,n,1),pdis(0:3,n,1)
     229             :     ENDDO
     230           0 :     WRITE(6,*) "Difference of interstitial:",pdis(0,0,1)
     231           0 :     IF (input%jspins.EQ.2) THEN
     232           0 :        WRITE (6,FMT=8010) iter,dis(2)
     233             : 8010   FORMAT (/,'----> distance of spin potential for it=',i3,':', f11.6,' mhtr/bohr**3')
     234           0 :        WRITE(6,*) "Details of potential differences for each atom type"
     235           0 :        WRITE(6,*) "Atom: total difference: difference of first sphhar"
     236           0 :        DO n=1,atoms%ntype
     237           0 :           WRITE(6,"(i5,' : ',f10.6,' : ',4f10.6)") n,pdis(4,n,2),pdis(0:3,n,2)
     238             :        ENDDO
     239           0 :        WRITE(6,*) "Difference of interstitial:",pdis(0,0,2)
     240             :     END IF
     241             :   CONTAINS
     242           0 :     SUBROUTINE priv_cdndif(m1,m2,m3,n1,n2,n3,s1,s2,film,rhsp,rhpw,rhv0,rhv1)
     243             :       IMPLICIT NONE
     244             :       INTEGER, INTENT (IN) :: m1,m2,m3,n1,n2,n3
     245             :       REAL,    INTENT (IN) :: s1,s2
     246             :       COMPLEX, INTENT (INOUT) :: rhpw(:,:,:)
     247             :       COMPLEX, INTENT (INOUT) :: rhv1(:,:,:,:,:)
     248             :       REAL,    INTENT (INOUT) :: rhsp(:,:,:,:,:)
     249             :       REAL,    INTENT (INOUT) :: rhv0(:,:,:,:)
     250             :       LOGICAL,INTENT (IN)     :: film
     251             : 
     252           0 :       rhsp(:,:,:,m1,n1)=s1*rhsp(:,:,:,m2,n2) +  s2*rhsp(:,:,:,m3,n3)
     253           0 :       rhpw(:,m1,n1)    =s1*rhpw(:,m2,n2)     +  s2*rhpw(:,m3,n3)
     254           0 :       IF (film) THEN
     255           0 :          rhv0(:,:,m1,n1)  =s1*rhv0(:,:,m2,n2)   + s2*rhv0(:,:,m3,n3)
     256           0 :          rhv1(:,:,:,m1,n1)=s1*rhv1(:,:,:,m2,n2) + s2*rhv1(:,:,:,m3,n3)
     257             :       ENDIF
     258           0 :     END SUBROUTINE priv_cdndif
     259             : 
     260             :   END SUBROUTINE potdis
     261             : END MODULE m_potdis

Generated by: LCOV version 1.13