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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_force_a4_add
       8             :   USE m_juDFT
       9             :   !     *****************************************************************
      10             :   !     calculates the force contribution from core-tails in addition
      11             :   !     to formula A4 of Yu, Singh & Krakauer minus a surface term that
      12             :   !     was included conveniently in force_a4.f
      13             :   !     It evaluates for the core density of each atom the term
      14             :   !     int [nabla rho_core]V d^3r in each region but the own muffin tin
      15             :   !     (i.e. the interstitial, the other muffin tins and the vacuum)
      16             :   !     (the corresponding formula is extended over the unit cell)
      17             :   !     Klueppelberg Sep'12
      18             :   !     *****************************************************************
      19             : 
      20             :   IMPLICIT NONE
      21             : 
      22             :   ! !     .. Local Array ..
      23             :   !       REAL   , ALLOCATABLE, PRIVATE, SAVE :: acoff(:,:),alpha(:,:)
      24             :   !       COMPLEX, ALLOCATABLE, PRIVATE, SAVE :: qpwcatom(:,:,:)
      25             :   !     .. Local Arrays ..
      26             :   REAL   , ALLOCATABLE, PUBLIC, SAVE :: force_a4_mt(:,:,:)
      27             :   COMPLEX, ALLOCATABLE, PUBLIC, SAVE :: force_a4_is(:,:,:)
      28             :   INTEGER, PUBLIC, SAVE :: f_level=3
      29             :   !     f_level == 0: Original force calculation
      30             :   !     f_level == 1: Forces from coretails calculated over whole unit cell
      31             :   !     f_level == 2: Kinetic energy surface term evaluated with IR functions
      32             :   !     f_level == 3: Surface term for density and potential discontinuity at the MT boundaries
      33             :   !     level 3 needs a large gmax cutoff, which backfires at the rest of the calculation
      34             :   !     TO DO: Read f_level from input file, right now, it is only set here
      35             : 
      36             : CONTAINS
      37             : 
      38             : 
      39             : 
      40           0 :   SUBROUTINE alloc_fa4_arrays(atoms,input)
      41             :     !     This subroutine allocates the arrays filled in cdnovlp.F
      42             :     !     so their content can be provided in totale.f, where the addition takes place
      43             :     USE m_types
      44             :     IMPLICIT NONE
      45             :     TYPE(t_input),INTENT(IN)   :: input
      46             :     TYPE(t_atoms),INTENT(IN)       :: atoms
      47             : 
      48             :     !     .. Scalar Arguments ..
      49             : 
      50           0 :     IF (.not.allocated(force_a4_mt)) THEN
      51           0 :        ALLOCATE ( force_a4_mt(3,atoms%ntype,input%jspins),force_a4_is(3,atoms%ntype,input%jspins) )
      52             :     END IF
      53             : 
      54           0 :   END SUBROUTINE alloc_fa4_arrays
      55             : 
      56             : 
      57             : 
      58           0 :   SUBROUTINE force_a4_add(&
      59             :        &                        atoms,input,&
      60           0 :        &                        force)
      61             :     !     This subroutine adds the coretail contribution of the forces to
      62             :     !     the atomic forces
      63             :     USE m_types
      64             :     IMPLICIT NONE
      65             :     TYPE(t_input),INTENT(IN)       :: input
      66             :     TYPE(t_atoms),INTENT(IN)       :: atoms
      67             :     REAL,INTENT(INOUT)             :: force(:,:,:)
      68             :     !     .. Local Scalars ..
      69             :     INTEGER :: jsp,n,dir
      70             :    
      71           0 :     CALL timestart("force_a4_add")
      72             : 
      73             :     ! The following film part needs to be implemented into cdnovlp, since
      74             :     ! the coretail force contribution has been relocated to cdnovlp.F
      75             :     !       CALL cpu_time(time1)
      76             :     ! !**********************************************************************
      77             :     ! !     in case of film calculations: calculate integral over vacuum
      78             :     ! !**********************************************************************
      79             :     !       force_a4_1d = czero
      80             :     !       force_a4_2d = czero
      81             :     !       IF (film.AND..not.odi%d1) THEN
      82             :     !       DO i = 1,max(nmzxy,nmz)
      83             :     !         z(i) = i*delz
      84             :     !       END DO ! i
      85             :     !       DO jsp = 1,jspins
      86             :     !       DO ivac = 1,nvac
      87             :     !       sumG = czero
      88             :     !       integrand2 = czero
      89             :     !       DO k = 2,nq3
      90             :     !         k1 = kv3(1,k)
      91             :     !         k2 = kv3(2,k)
      92             :     !         k3 = kv3(3,k)
      93             :     !         ind2 = ig2(k)
      94             :     ! !         IF (ind2.EQ.1) CYCLE
      95             :     ! 
      96             :     !         CALL spgrot(
      97             :     !      >              nop,symor,tpi,mrot,tau,invtab,
      98             :     !      >              kv3(:,k),
      99             :     !      <              kr3d,phas3d)
     100             :     ! 
     101             :     !         CALL spgrot(
     102             :     !      >              nop2,symor,tpi,mrot,tau,invtab,
     103             :     !      >              (/k1,k2,0/),
     104             :     !      <              kr2d,phas2d)
     105             :     ! 
     106             :     !         carg = ci * k3 * bmat(3,3) ! k3 in external coordinates
     107             :     ! 
     108             :     !         pot = 0.0
     109             :     !         IF (ind2.EQ.1) THEN
     110             :     !           pot(1:nmz  ) = vz( 1:nmz         ,ivac,jsp)
     111             :     !         ELSE
     112             :     !           pot(1:nmzxy) = vxy(1:nmzxy,ind2-1,ivac,jsp)
     113             :     !         END IF
     114             :     ! 
     115             :     ! 
     116             :     !         DO j = 1,nop
     117             :     !           kp = (k-1)*nop + j
     118             :     !         DO jp = 1,nop2
     119             :     !           IF (kr3d(1,j).NE.kr2d(1,jp)) CYCLE
     120             :     !           IF (kr3d(2,j).NE.kr2d(2,jp)) CYCLE
     121             :     !           sumG(:,kp) = sumG(:,kp) +exp(carg*z(:))*pot(:)*phas2d(jp)/nop2
     122             :     !         END DO ! jp 2d operations
     123             :     !         END DO ! j operations
     124             :     ! 
     125             :     !         IF (k.EQ.nq3) THEN
     126             :     !         DO n = 1,ntype
     127             :     !           DO j = 1,nop
     128             :     !             kp = (k-1)*nop + j
     129             :     !             DO dir = 1,3
     130             :     !               integrand2(dir,1:mshd,n) = integrand2(dir,1:mshd,n)
     131             :     !      +                                 + dqpwcatom(dir,kp,n,jsp)
     132             :     !      *                                 * sumG(1:mshd,kp)
     133             :     !             END DO ! dir ections
     134             :     !           END DO ! j operations
     135             :     !         END DO ! n types of atoms
     136             :     !         END IF
     137             :     ! 
     138             :     !       END DO ! k stars
     139             :     ! 
     140             :     !       DO n = 1,ntype
     141             :     !       DO dir = 1,3
     142             :     !        CALL qsf(delz, real(integrand2(dir,:,n)),integral2(1,dir),mshd,0)
     143             :     !        CALL qsf(delz,aimag(integrand2(dir,:,n)),integral2(2,dir),mshd,0)
     144             :     !         integral2(:,dir) = integral2(:,dir) * area
     145             :     !         force_a4_2d(dir,n,jsp) = force_a4_2d(dir,n,jsp)
     146             :     !      +                         + integral2(1,dir) + ci* integral2(2,dir)
     147             :     !       END DO ! dir ections
     148             :     !       END DO ! n types of atoms
     149             :     ! 
     150             :     !       END DO ! ivac vacuum regions (1 or 2)
     151             :     !       END DO ! jsp spins
     152             :     !       END IF
     153             :     !       CALL cpu_time(time2)
     154             :     !       WRITE (*,*) 'time to calculate lower dimensional contributions:',
     155             :     !      + time2-time1
     156             : 
     157             :     !**********************************************************************
     158             :     !     add results to total force
     159             :     !**********************************************************************
     160           0 :     DO jsp = 1,input%jspins
     161           0 :        DO n = 1,atoms%ntype
     162             : 
     163           0 :           IF (atoms%l_geo(n)) THEN
     164             : 
     165             :              force(:,n,jsp) = force(:,n,jsp) + real(force_a4_is(:,n,jsp))&
     166           0 :                   &                                  + real(force_a4_mt(:,n,jsp))
     167             :              !      +                                  + real(force_a4_2d(:,n,jsp)) ! is calculated only if film.and..not.odi%d1, otherwise 0
     168             :              !      +                                  + real(force_a4_1d(:,n,jsp)) ! is calculated only if odi%d1, otherwise 0
     169             : 
     170             :              !       write to out-file
     171           0 :              WRITE (6,FMT=8010) n
     172           0 :              WRITE (6,FMT=8020) ((force_a4_is(dir,n,jsp)),dir=1,3) ! 8020
     173           0 :              WRITE (6,FMT=8015) n
     174           0 :              WRITE (6,FMT=8020) ((force_a4_mt(dir,n,jsp)),dir=1,3) ! 8020
     175             : 8010         FORMAT (' FORCES: IS ADDITION TO EQUATION A4 FOR ATOM TYPE',i4)
     176             : 8015         FORMAT (' FORCES: MT ADDITION TO EQUATION A4 FOR ATOM TYPE',i4)
     177             :              !  8025   FORMAT (' FORCES: VACUUM ADD. TO EQUATION A4 FOR ATOM TYPE',i4)
     178             :              !  8020   FORMAT (' FX_A4=',2f19.15,' FY_A4=',2f19.15,' FZ_A4=',2f19.15)
     179             : 8020         FORMAT (' FX_A4=',2f10.6,' FY_A4=',2f10.6,' FZ_A4=',2f10.6)
     180             : 8070         FORMAT (' FX_A4=',2f20.14,' FY_A4=',2f20.14,' FZ_A4=',2f20.14)
     181             : 
     182             :           END IF ! atoms%l_geo(n)
     183             : 
     184             :        END DO ! n types of atoms
     185             :     END DO ! jsp spins
     186             :     
     187           0 :     DEALLOCATE ( force_a4_mt,force_a4_is )
     188           0 :     call timestop("force_a4_add")
     189             : 
     190           0 :   END SUBROUTINE force_a4_add
     191             : 
     192             : END MODULE m_force_a4_add

Generated by: LCOV version 1.13