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 : USE m_types
10 : !------------------------------------------------------------------------------
11 : ! Adds the force contribution from core-tails in addition to formula A4 of Yu,
12 : ! Singh & Krakauer minus a surface term that was included conveniently in
13 : ! force_a4.f90. The contributions are calculated in cdnovlp.F90.
14 : !
15 : ! It evaluates the term int [nabla rho_core]V d^3r for the core density of each
16 : ! atom in each region but its own Muffin Tin (i.e. the interstitial and the
17 : ! other MT).
18 : ! The corresponding formula is extended over the unit cell.
19 : ! Klueppelberg Sep'12 (force level 1)
20 : !------------------------------------------------------------------------------
21 :
22 : IMPLICIT NONE
23 :
24 : REAL , ALLOCATABLE, PUBLIC, SAVE :: force_a4_mt(:,:,:)
25 : COMPLEX, ALLOCATABLE, PUBLIC, SAVE :: force_a4_is(:,:,:)
26 :
27 : CONTAINS
28 2 : SUBROUTINE alloc_fa4_arrays(atoms,input)
29 : ! This subroutine allocates the arrays filled in cdnovlp.F90
30 : ! so their content can be provided in totale.f90, where the addition takes place
31 :
32 : TYPE(t_input), INTENT(IN) :: input
33 : TYPE(t_atoms), INTENT(IN) :: atoms
34 :
35 2 : CALL timestart("force_a4_alloc")
36 :
37 0 : ALLOCATE(force_a4_mt(3,atoms%ntype,input%jspins), &
38 14 : force_a4_is(3,atoms%ntype,input%jspins))
39 :
40 2 : CALL timestop("force_a4_alloc")
41 :
42 2 : END SUBROUTINE alloc_fa4_arrays
43 :
44 1 : SUBROUTINE force_a4_add(atoms,input,results)
45 : ! This subroutine adds the coretail contribution to the atomic forces.
46 :
47 : USE m_constants
48 :
49 : TYPE(t_input), INTENT(IN) :: input
50 : TYPE(t_atoms), INTENT(IN) :: atoms
51 : TYPE(t_results), INTENT(INOUT) :: results
52 :
53 : INTEGER :: jsp,n,dir
54 :
55 1 : COMPLEX :: force_a4_mt_loc(3,atoms%ntype,input%jspins)
56 :
57 1 : CALL timestart("force_a4_add")
58 :
59 : ! The following film part needs to be implemented into cdnovlp, since
60 : ! the coretail force contribution has been relocated to cdnovlp.F
61 : ! CALL cpu_time(time1) TODO: Readd this film part.
62 : ! !**********************************************************************
63 : ! ! in case of film calculations: calculate integral over vacuum
64 : ! !**********************************************************************
65 : ! force_a4_1d = czero
66 : ! force_a4_2d = czero
67 : ! IF (film.AND..not.odi%d1) THEN
68 : ! DO i = 1,max(nmzxy,nmz)
69 : ! z(i) = i*delz
70 : ! END DO ! i
71 : ! DO jsp = 1,jspins
72 : ! DO ivac = 1,nvac
73 : ! sumG = czero
74 : ! integrand2 = czero
75 : ! DO k = 2,nq3
76 : ! k1 = kv3(1,k)
77 : ! k2 = kv3(2,k)
78 : ! k3 = kv3(3,k)
79 : ! ind2 = ig2(k)
80 : ! ! IF (ind2.EQ.1) CYCLE
81 : !
82 : ! CALL spgrot(
83 : ! > nop,symor,tpi,mrot,tau,invtab,
84 : ! > kv3(:,k),
85 : ! < kr3d,phas3d)
86 : !
87 : ! CALL spgrot(
88 : ! > nop2,symor,tpi,mrot,tau,invtab,
89 : ! > (/k1,k2,0/),
90 : ! < kr2d,phas2d)
91 : !
92 : ! carg = ci * k3 * bmat(3,3) ! k3 in external coordinates
93 : !
94 : ! pot = 0.0
95 : ! IF (ind2.EQ.1) THEN
96 : ! pot(1:nmz ) = vz( 1:nmz ,ivac,jsp)
97 : ! ELSE
98 : ! pot(1:nmzxy) = vxy(1:nmzxy,ind2-1,ivac,jsp)
99 : ! END IF
100 : !
101 : !
102 : ! DO j = 1,nop
103 : ! kp = (k-1)*nop + j
104 : ! DO jp = 1,nop2
105 : ! IF (kr3d(1,j).NE.kr2d(1,jp)) CYCLE
106 : ! IF (kr3d(2,j).NE.kr2d(2,jp)) CYCLE
107 : ! sumG(:,kp) = sumG(:,kp) +exp(carg*z(:))*pot(:)*phas2d(jp)/nop2
108 : ! END DO ! jp 2d operations
109 : ! END DO ! j operations
110 : !
111 : ! IF (k.EQ.nq3) THEN
112 : ! DO n = 1,ntype
113 : ! DO j = 1,nop
114 : ! kp = (k-1)*nop + j
115 : ! DO dir = 1,3
116 : ! integrand2(dir,1:mshd,n) = integrand2(dir,1:mshd,n)
117 : ! + + dqpwcatom(dir,kp,n,jsp)
118 : ! * * sumG(1:mshd,kp)
119 : ! END DO ! directions
120 : ! END DO ! j operations
121 : ! END DO ! n types of atoms
122 : ! END IF
123 : !
124 : ! END DO ! k stars
125 : !
126 : ! DO n = 1,ntype
127 : ! DO dir = 1,3
128 : ! CALL qsf(delz, real(integrand2(dir,:,n)),integral2(1,dir),mshd,0)
129 : ! CALL qsf(delz,aimag(integrand2(dir,:,n)),integral2(2,dir),mshd,0)
130 : ! integral2(:,dir) = integral2(:,dir) * area
131 : ! force_a4_2d(dir,n,jsp) = force_a4_2d(dir,n,jsp)
132 : ! + + integral2(1,dir) + ci* integral2(2,dir)
133 : ! END DO ! dir ections
134 : ! END DO ! n types of atoms
135 : !
136 : ! END DO ! ivac vacuum regions (1 or 2)
137 : ! END DO ! jsp spins
138 : ! END IF
139 : ! CALL cpu_time(time2)
140 : ! WRITE (*,*) 'time to calculate lower dimensional contributions:',
141 : ! + time2-time1
142 :
143 : ! Add the results to the total force.
144 :
145 14 : force_a4_mt_loc=force_a4_mt*CMPLX(1.0,0.0)
146 2 : DO jsp = 1, input%jspins
147 5 : DO n = 1, atoms%ntype
148 4 : IF (atoms%l_geo(n)) THEN
149 :
150 12 : results%force(:,n,jsp) = results%force(:,n,jsp) + real(force_a4_is(:,n,jsp)) + real(force_a4_mt(:,n,jsp))
151 : ! + real(force_a4_2d(:,n,jsp)) ! is calculated only if film.and..not.odi%d1, otherwise 0
152 : ! + real(force_a4_1d(:,n,jsp)) ! is calculated only if odi%d1, otherwise 0
153 :
154 : ! Write force contributions to out file.
155 3 : WRITE (oUnit,FMT=8010) n
156 12 : WRITE (oUnit,FMT=8020) ((force_a4_is(dir,n,jsp)),dir=1,3) ! 8020
157 3 : WRITE (oUnit,FMT=8015) n
158 12 : WRITE (oUnit,FMT=8070) ((force_a4_mt_loc(dir,n,jsp)),dir=1,3) ! 8070
159 : 8010 FORMAT (' FORCES: IS ADDITION TO EQUATION A4 FOR ATOM TYPE',i4)
160 : 8015 FORMAT (' FORCES: MT ADDITION TO EQUATION A4 FOR ATOM TYPE',i4)
161 : !8025 FORMAT (' FORCES: VACUUM ADD. TO EQUATION A4 FOR ATOM TYPE',i4)
162 : !8020 FORMAT (' FX_A4=',2f19.15,' FY_A4=',2f19.15,' FZ_A4=',2f19.15)
163 : 8020 FORMAT (' FX_IS=',2f10.6,' FY_IS=',2f10.6,' FZ_IS=',2f10.6)
164 : 8070 FORMAT (' FX_MT=',2f10.6,' FY_MT=',2f10.6,' FZ_MT=',2f10.6)
165 :
166 : END IF ! atoms%l_geo(n)
167 :
168 : END DO ! n types of atoms
169 : END DO ! jsp spins
170 :
171 1 : DEALLOCATE ( force_a4_mt, force_a4_is )
172 1 : CALL timestop("force_a4_add")
173 :
174 1 : END SUBROUTINE force_a4_add
175 :
176 : END MODULE m_force_a4_add
|
