Branch data Line data Source code
1 : : !!****m* ABINIT/m_dens
2 : : !! NAME
3 : : !! m_dens
4 : : !!
5 : : !! FUNCTION
6 : : !! Module containing the definition of the constrained_dft_t data type and methods used to handle it,
7 : : !! and also includes the computation of integrated atomic charge and magnetization, as well as Hirshfeld charges.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 1998-2026 ABINIT group (MT,ILuk,MVer,EB,SPr)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : MODULE m_dens
24 : :
25 : : use defs_basis
26 : : use m_errors
27 : : use m_abicore
28 : : use m_xmpi
29 : : use m_splines
30 : :
31 : : use defs_abitypes, only : MPI_type
32 : : use m_fft, only : fourdp,fftpac
33 : : use m_time, only : timab
34 : : use m_numeric_tools, only : wrap2_zero_one, geteuler
35 : : use m_io_tools, only : open_file
36 : : use m_geometry, only : dist2, xcart2xred, metric, vcart2ylm, cart2spinaxis
37 : : use m_mpinfo, only : ptabs_fourdp
38 : : use m_atomdata
39 : : use m_dtset
40 : :
41 : : implicit none
42 : :
43 : : private
44 : :
45 : : public :: dens_hirsh ! Compute the Hirshfeld charges
46 : : public :: add_atomic_fcts ! Add atomic functions to real space function
47 : : public :: constrained_dft_ini ! Initialize the constrained_dft datastructure
48 : : public :: constrained_dft_free ! Free the constrained_dft datastructure
49 : : public :: constrained_residual ! Recompute the potential residual, to account for constraints
50 : : public :: mag_penalty ! Compute the potential corresponding to constrained magnetic moments (using add_atomic_fcts) with the penalty function.
51 : : public :: mag_penalty_e ! Compute the energy corresponding to constrained magnetic moments.
52 : : public :: calcdenmagsph ! Compute integral of total density and magnetization inside spheres around atoms.
53 : : public :: prtdenmagsph ! Print integral of total density and magnetization inside spheres around atoms.
54 : : public :: magmom_to_d2 ! Integrates the magnetic moments (total & local ones) into the ddb array.
55 : : public :: fatsph_recip ! Compute atom centered spheres in reciprocal space
56 : : public :: calmaxdifmag ! Compute the maximum magnetization and the maximum change in magnetization between two steps.
57 : : !!***
58 : :
59 : : !----------------------------------------------------------------------
60 : :
61 : : !!****t* m_dens/constrained_dft_t
62 : : !! NAME
63 : : !! constrained_dft_t
64 : : !!
65 : : !! FUNCTION
66 : : !! Structure gathering the relevant information for constrained DFT calculations
67 : : !!
68 : : !! SOURCE
69 : :
70 : : type,public :: constrained_dft_t
71 : :
72 : : !scalars
73 : : integer :: natom ! Number of atoms
74 : : integer :: nfftf ! Number of FFT grid points (for this processor) for the "fine" grid
75 : : integer :: nspden ! Number of spin-density components
76 : : integer :: ntypat ! Number of type of atoms
77 : :
78 : : real(dp) :: magcon_lambda ! Strength of the atomic spherical constraint
79 : : real(dp) :: ratsm ! Smearing width for ratsph
80 : :
81 : : integer :: magconon ! Turn on the penalty function constraint instead of the more powerful constrainedDFT algorithm
82 : :
83 : : !arrays
84 : :
85 : : integer :: ngfftf(18) ! Number of FFT grid points (for this processor) for the "fine" grid
86 : :
87 : : integer,allocatable :: typat(:)
88 : : ! typat(natom)
89 : : ! Type of each natom
90 : :
91 : : integer,allocatable :: constraint_kind(:)
92 : : ! constraint_kind(ntypat)
93 : : ! Constraint kind to be applied to each type of atom. See corresponding input variable
94 : :
95 : : real(dp) :: rprimd(3,3)
96 : : ! Direct lattice vectors, Bohr units.
97 : :
98 : : real(dp),allocatable :: chrgat(:)
99 : : ! chrgat(natom)
100 : : ! Target charge for each atom. Not always used, it depends on the value of constraint_kind
101 : :
102 : : real(dp),allocatable :: intgf2(:,:)
103 : : ! intgf2(natom,natom)
104 : : ! Overlap of the spherical integrating functions, for each atom.
105 : : ! Initialized using some xred values, will not change during the SCF cycles, except for exotic algorithms, not in production,
106 : :
107 : : real(dp),allocatable :: ratsph(:)
108 : : ! ratsph(ntypat)
109 : : ! Radius of the atomic sphere for each type of atom
110 : :
111 : : real(dp),allocatable :: spinat(:,:)
112 : : ! spinat(3,natom)
113 : : ! Target magnetization for each atom. Possibly only the direction or the magnitude, depending on constraint_kind
114 : :
115 : : real(dp),allocatable :: ziontypat(:)
116 : : ! ziontypat(ntypat)
117 : : ! Ionic charge, per type of atom
118 : :
119 : : real(dp),allocatable :: znucl(:)
120 : : ! znucl(ntypat)
121 : : !
122 : :
123 : : real(dp) :: spinaxis(3)
124 : : ! spinaxis(3)
125 : : ! Spin quantization axis
126 : :
127 : : end type constrained_dft_t
128 : :
129 : : !!***
130 : :
131 : : CONTAINS
132 : :
133 : : !----------------------------------------------------------------------
134 : :
135 : : !!****f* m_dens/dens_hirsh
136 : : !! NAME
137 : : !! dens_hirsh
138 : : !!
139 : : !! FUNCTION
140 : : !! Compute the Hirshfeld charges
141 : : !!
142 : : !! INPUTS
143 : : !! mpoint=Maximum number of points in radial meshes.
144 : : !! radii(mpoint, ntypat)=Radial meshes for each type
145 : : !! aeden(mpoint, nytpat)=All-electron densities.
146 : : !! npoint(ntypat)=The number of the last point with significant density is stored in npoint(itypat)
147 : : !! minimal_den=Tolerance on the minum value of the density
148 : : !! grid_den(nrx,nry,nrz)= density on the grid
149 : : !! natom = number of atoms in the unit cell
150 : : !! nrx,nry,nrz= number of points in the grid for the three directions
151 : : !! ntypat=number of types of atoms in unit cell.
152 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
153 : : !! typat(natom)=type of each atom
154 : : !! xcart(3,natom) = different positions of the atoms in the unit cell
155 : : !! zion=(ntypat)gives the ionic charge for each type of atom
156 : : !! prtcharge=1 to write the Hirshfeld charge decomposition
157 : : !!
158 : : !! OUTPUT
159 : : !! hcharge(natom), hden(natom), hweight(natom)= Hirshfeld charges, densities, weights.
160 : : !!
161 : : !! SOURCE
162 : :
163 : 1 : subroutine dens_hirsh(mpoint,radii,aeden,npoint,minimal_den,grid_den, &
164 : 1 : natom,nrx,nry,nrz,ntypat,rprimd,xcart,typat,zion,prtcharge,hcharge,hden,hweight)
165 : :
166 : : !Arguments ------------------------------------
167 : : !scalars
168 : : integer,intent(in) :: natom,nrx,nry,nrz,ntypat,prtcharge,mpoint
169 : : real(dp),intent(in) :: minimal_den
170 : : !arrays
171 : : integer,intent(in) :: typat(natom),npoint(ntypat)
172 : : real(dp),intent(in) :: grid_den(nrx,nry,nrz),rprimd(3,3),zion(ntypat)
173 : : real(dp),intent(in) :: xcart(3,natom)
174 : : real(dp),intent(in) :: radii(mpoint,ntypat),aeden(mpoint,ntypat)
175 : : real(dp),intent(out) :: hcharge(natom),hden(natom),hweight(natom)
176 : :
177 : : !Local variables -------------------------
178 : : !scalars
179 : : integer :: i1,i2,i3,iatom,icell,igrid,ii,inmax,inmin,istep,itypat
180 : : integer :: k1,k2,k3,mcells,nfftot,ngoodpoints,npt
181 : : real(dp) :: aa,bb,coeff1,coeff2,coeff3,den,factor,h_inv,hh,maxrad
182 : : real(dp) :: rr,rr2,total_charge,total_weight,total_zion,ucvol
183 : : real(dp) :: yp1,ypn
184 : : !arrays
185 : : integer :: highest(3),lowest(3)
186 : 1 : integer,allocatable :: ncells(:)
187 : : real(dp) :: coordat(3),coord23_1,coord23_2,coord23_3,diff1,diff2,diff3,gmet(3,3),gprimd(3,3),rmet(3,3)
188 : : real(dp) :: vperp(3),width(3)
189 : 1 : real(dp),allocatable :: coord1(:,:),local_den(:,:,:,:)
190 : 1 : real(dp),allocatable :: step(:,:),sum_den(:,:,:)
191 : 1 : real(dp),allocatable :: xcartcells(:,:,:),xred(:,:),yder2(:)
192 : : ! *********************************************************************
193 : :
194 : : !1. Read the 1D all-electron atomic files
195 : : !Store the radii in radii(:,itypat), and the all-electron
196 : : !densities in aeden(:,itypat). The number of the last
197 : : !point with significant density is stored in npoint(itypat)
198 : :
199 : : !2. Compute the list of atoms that are sufficiently close to the cell
200 : :
201 : 1 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
202 : 1 : nfftot=nrx*nry*nrz
203 : :
204 [ - + + - : 3 : ABI_MALLOC(xred,(3,natom))
+ - + - -
+ - + ]
205 : 1 : call xcart2xred(natom,rprimd,xcart,xred)
206 : :
207 : : !Compute the widths of the cell
208 : : !First width : perpendicular vector length
209 [ + + ]: 4 : vperp(:)=rprimd(:,1)-rprimd(:,2)*rmet(1,2)/rmet(2,2) -rprimd(:,3)*rmet(1,3)/rmet(3,3)
210 [ + + ]: 4 : width(1)=sqrt(dot_product(vperp,vperp))
211 : : !Second width
212 [ + + ]: 4 : vperp(:)=rprimd(:,2)-rprimd(:,1)*rmet(2,1)/rmet(1,1) -rprimd(:,3)*rmet(2,3)/rmet(3,3)
213 [ + + ]: 4 : width(2)=sqrt(dot_product(vperp,vperp))
214 : : !Third width
215 [ + + ]: 4 : vperp(:)=rprimd(:,3)-rprimd(:,1)*rmet(3,1)/rmet(1,1) -rprimd(:,2)*rmet(3,2)/rmet(2,2)
216 [ + + ]: 4 : width(3)=sqrt(dot_product(vperp,vperp))
217 : :
218 : : !Compute the number of cells that will make up the supercell
219 [ - + + - : 3 : ABI_MALLOC(ncells,(natom))
+ - + - -
+ - + ]
220 : 1 : mcells=0
221 [ + + ]: 3 : do iatom=1,natom
222 : 2 : itypat=typat(iatom)
223 : 2 : maxrad=radii(npoint(itypat),itypat)
224 : : ! Compute the lower and higher indices of the supercell
225 : : ! for this atom
226 [ + + ]: 8 : do ii=1,3
227 [ + - ]: 6 : lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii))
228 [ + - ]: 8 : highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1)
229 : : ! Next coding, still incorrect
230 : : ! lowest(ii)=floor(xred(ii,iatom)-maxrad/width(ii))-1
231 : : ! highest(ii)=ceiling(xred(ii,iatom)+maxrad/width(ii))+1
232 : : ! Old coding incorrect
233 : : ! lowest(ii)=ceiling(-xred(ii,iatom)-maxrad/width(ii))
234 : : ! highest(ii)=floor(-xred(ii,iatom)+maxrad/width(ii)+1)
235 : : end do
236 : : ncells(iatom)=(highest(1)-lowest(1)+1)* &
237 : : & (highest(2)-lowest(2)+1)* &
238 : 3 : & (highest(3)-lowest(3)+1)
239 : : ! DEBUG
240 : : ! write(std_out,*)' maxrad=',maxrad
241 : : ! write(std_out,*)' lowest(:)=',lowest(:)
242 : : ! write(std_out,*)' highest(:)=',highest(:)
243 : : ! write(std_out,*)' ncells(iatom)=',ncells(iatom)
244 : : ! ENDDEBUG
245 : : end do
246 [ + + ]: 3 : mcells=maxval(ncells(:))
247 : :
248 : : !Compute, for each atom, the set of image atoms in the whole supercell
249 [ - + + - : 4 : ABI_MALLOC(xcartcells,(3,mcells,natom))
- + + - +
- + - - +
- + ]
250 [ + + ]: 3 : do iatom=1,natom
251 : 2 : itypat=typat(iatom)
252 : 2 : maxrad=radii(npoint(itypat),itypat)
253 : : ! Compute the lower and higher indices of the supercell
254 : : ! for this atom
255 : :
256 [ + + ]: 8 : do ii=1,3
257 [ + - ]: 6 : lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii))
258 [ + - ]: 8 : highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1)
259 : : end do
260 : 2 : icell=0
261 [ + + ]: 16 : do i1=lowest(1),highest(1)
262 [ + + ]: 104 : do i2=lowest(2),highest(2)
263 [ + + ]: 739 : do i3=lowest(3),highest(3)
264 : 637 : icell=icell+1
265 [ + + ]: 2637 : xcartcells(:,icell,iatom)=xcart(:,iatom)+i1*rprimd(:,1)+i2*rprimd(:,2)+i3*rprimd(:,3)
266 : : end do
267 : : end do
268 : : end do
269 : : end do
270 : :
271 : : !Compute, for each atom, the all-electron pro-atom density
272 : : !at each point in the primitive cell
273 [ - + + - : 6 : ABI_MALLOC(local_den,(nrx,nry,nrz,natom))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
274 [ - + + - : 3 : ABI_MALLOC(step,(2,mpoint))
+ - + - -
+ - + ]
275 [ - + + - : 3 : ABI_MALLOC(yder2,(mpoint))
+ - + - -
+ - + ]
276 [ - + + - : 3 : ABI_MALLOC(coord1,(3,nrx))
+ - + - -
+ - + ]
277 : 1 : coeff1=one/nrx
278 : 1 : coeff2=one/nry
279 : 1 : coeff3=one/nrz
280 : :
281 [ + + ]: 3 : do iatom=1,natom
282 : 2 : itypat=typat(iatom)
283 : 2 : npt=npoint(itypat)
284 : 2 : maxrad=radii(npt,itypat)
285 : : ! write(std_out,*)
286 : : ! write(std_out,'(a,i4)' )' hirsh : accumulating density for atom ',iatom
287 : : ! write(std_out,*)' ncells(iatom)=',ncells(iatom)
288 [ + + ]: 829 : do istep=1,npt-1
289 : 827 : step(1,istep)=radii(istep+1,itypat) - radii(istep,itypat)
290 : 829 : step(2,istep)=one/step(1,istep)
291 : : end do
292 : : ! Approximate first derivative for small radii
293 : 2 : yp1=(aeden(2,itypat)-aeden(1,itypat))/(radii(2,itypat)-radii(1,itypat))
294 : 2 : ypn=zero
295 : 2 : call spline(radii(1:npt,itypat),aeden(1:npt,itypat),npt,yp1,ypn,yder2)
296 : :
297 [ + + + + : 95978 : local_den(:,:,:,iatom)=zero
+ + ]
298 : :
299 : : ! Big loop on the cells
300 [ + + ]: 640 : do icell=1,ncells(iatom)
301 : : ! write(std_out,*)' icell=',icell
302 [ + + ]: 2548 : coordat(:)=xcartcells(:,icell,iatom)
303 : :
304 : : ! Big loop on the grid points
305 [ + + ]: 23569 : do k1 = 1,nrx
306 [ + + ]: 92365 : coord1(:,k1)=rprimd(:,1)*(k1-1)*coeff1
307 : : end do
308 [ + + ]: 23571 : do k3 = 1, nrz
309 [ + + ]: 849121 : do k2 = 1, nry
310 : 825552 : coord23_1=rprimd(1,2)*(k2-1)*coeff2+rprimd(1,3)*(k3-1)*coeff3-coordat(1)
311 : 825552 : coord23_2=rprimd(2,2)*(k2-1)*coeff2+rprimd(2,3)*(k3-1)*coeff3-coordat(2)
312 : 825552 : coord23_3=rprimd(3,2)*(k2-1)*coeff2+rprimd(3,3)*(k3-1)*coeff3-coordat(3)
313 [ + + ]: 30568356 : do k1 = 1, nrx
314 : 29719872 : diff1=coord1(1,k1)+coord23_1
315 : 29719872 : diff2=coord1(2,k1)+coord23_2
316 : 29719872 : diff3=coord1(3,k1)+coord23_3
317 : 29719872 : rr2=diff1**2+diff2**2+diff3**2
318 [ + + ]: 30545424 : if(rr2<maxrad**2)then
319 : :
320 : 2266812 : rr=sqrt(rr2)
321 : : ! Find the index of the radius by bissection
322 [ + + ]: 2266812 : if (rr < radii(1,itypat)) then
323 : : ! Linear extrapolation
324 : : den=aeden(1,itypat)+(rr-radii(1,itypat))/(radii(2,itypat)-radii(1,itypat))&
325 : 2 : & *(aeden(2,itypat)-aeden(1,itypat))
326 : : else
327 : : ! Use the spline interpolation
328 : : ! Find the index of the radius by bissection
329 : : inmin=1
330 : : inmax=npt
331 : 22255544 : igrid=1
332 : : do
333 [ + + ]: 22255544 : if(inmax-inmin==1)exit
334 : 19988734 : igrid=(inmin+inmax)/2
335 [ + + ]: 22255544 : if(rr>=radii(igrid,itypat))then
336 : : inmin=igrid
337 : : else
338 : 6662810 : inmax=igrid
339 : : end if
340 : : end do
341 : 2266810 : igrid=inmin
342 : : ! write(std_out,*)' igrid',igrid
343 : :
344 : 2266810 : hh=step(1,igrid)
345 : 2266810 : h_inv=step(2,igrid)
346 : 2266810 : aa= (radii(igrid+1,itypat)-rr)*h_inv
347 : 2266810 : bb= (rr-radii(igrid,itypat))*h_inv
348 : : den = aa*aeden(igrid,itypat) + bb*aeden(igrid+1,itypat) &
349 : : & +( (aa*aa*aa-aa)*yder2(igrid) &
350 : 2266810 : & +(bb*bb*bb-bb)*yder2(igrid+1) ) *hh*hh*sixth
351 : : end if ! Select small radius or spline
352 : :
353 : 2266812 : local_den(k1,k2,k3,iatom)=local_den(k1,k2,k3,iatom)+den
354 : : end if ! dist2<maxrad
355 : :
356 : : end do ! k1
357 : : end do ! k2
358 : : end do ! k3
359 : :
360 : : end do ! icell
361 : : end do ! iatom
362 : :
363 : : !Compute, the total all-electron density at each point in the primitive cell
364 [ - + + - : 5 : ABI_MALLOC(sum_den,(nrx,nry,nrz))
- + + - -
+ + - + -
+ - - + -
+ ]
365 [ + + + + : 47989 : sum_den(:,:,:)=zero
+ + ]
366 [ + + ]: 3 : do iatom=1,natom
367 [ + + + + : 95979 : sum_den(:,:,:)=sum_den(:,:,:)+local_den(:,:,:,iatom)
+ + ]
368 : : end do
369 : :
370 : : !Accumulate the integral of the density, to get Hirshfeld charges
371 : : !There is a minus sign because the electron has a negative charge
372 : : ngoodpoints = 0
373 [ + + ]: 3 : hcharge(:)=zero
374 [ + + ]: 3 : hweight(:)=zero
375 [ + + ]: 37 : do k3=1,nrz
376 [ + + ]: 1333 : do k2=1,nry
377 [ + + ]: 47988 : do k1=1,nrx
378 : : ! Use minimal_den in order to avoid divide by zero
379 [ + - ]: 47952 : if (abs(sum_den(k1,k2,k3)) > minimal_den) then
380 : : ngoodpoints = ngoodpoints+1
381 : 46656 : factor=grid_den(k1,k2,k3)/(sum_den(k1,k2,k3)+minimal_den)
382 [ + + ]: 139968 : do iatom=1,natom
383 : 93312 : hden(iatom)=hden(iatom)+local_den(k1,k2,k3,iatom)
384 : 93312 : hcharge(iatom)=hcharge(iatom)-local_den(k1,k2,k3,iatom)*factor
385 : 139968 : hweight(iatom)=hweight(iatom)+local_den(k1,k2,k3,iatom)/(sum_den(k1,k2,k3)+minimal_den)
386 : : end do
387 : : end if
388 : : end do
389 : : end do
390 : : end do
391 : :
392 : : !DEBUG
393 : : !do iatom=1,natom
394 : : !write(std_out,'(i9,3es17.6)' )iatom,hden(iatom),hcharge(iatom),hweight(iatom)
395 : : !end do
396 : : !ENDDEBUG
397 : :
398 [ + + ]: 3 : hcharge(:)=hcharge(:)*ucvol/dble(nfftot)
399 [ + + ]: 3 : hweight(:)=hweight(:)/dble(nfftot)
400 : :
401 : : !Check on the total charge
402 [ + + ]: 3 : total_zion=sum(zion(typat(1:natom)))
403 [ + + ]: 3 : total_charge=sum(hcharge(1:natom))
404 : : total_weight=sum(hweight(1:natom))
405 : :
406 : : !DEBUG
407 : : !write(std_out,*)' ngoodpoints = ', ngoodpoints, ' out of ', nfftot
408 : : !write(std_out,*)' total_weight=',total_weight
409 : : !write(std_out,*)' total_weight=',total_weight
410 : : !ENDDEBUG
411 : :
412 : : !Output
413 [ + - ]: 1 : if (prtcharge == 1) then
414 : 1 : write(std_out,*)
415 : 1 : write(std_out,*)' Hirshfeld analysis'
416 : 1 : write(std_out,*)' Atom Zion Electron Charge Net charge '
417 : 1 : write(std_out,*)
418 [ + + ]: 3 : do iatom=1,natom
419 : : write(std_out,'(i9,3es17.6)' )&
420 : 3 : & iatom,zion(typat(iatom)),hcharge(iatom),hcharge(iatom)+zion(typat(iatom))
421 : : end do
422 : 1 : write(std_out,*)
423 : 1 : write(std_out,'(a,3es17.6)')' Total',total_zion,total_charge,total_charge+total_zion
424 : 1 : write(std_out,*)
425 : : end if
426 : :
427 [ - + ]: 1 : ABI_FREE(coord1)
428 [ - + ]: 1 : ABI_FREE(local_den)
429 [ - + ]: 1 : ABI_FREE(ncells)
430 [ - + ]: 1 : ABI_FREE(step)
431 [ - + ]: 1 : ABI_FREE(sum_den)
432 [ - + ]: 1 : ABI_FREE(xcartcells)
433 [ - + ]: 1 : ABI_FREE(xred)
434 [ - + ]: 1 : ABI_FREE(yder2)
435 : :
436 : 1 : end subroutine dens_hirsh
437 : : !!***
438 : :
439 : : !!****f* m_dens/add_atomic_fcts
440 : : !! NAME
441 : : !! add_atomic_fcts
442 : : !!
443 : : !! FUNCTION
444 : : !! This routine is called to assemble the atomic spherical functions, and, if option/=0,
445 : : !! to add it to some input function (usually an input potential residual).
446 : : !! The contributions from each atomic sphere are governed by parameters coeff_constr_dft, input to the present routine.
447 : : !!
448 : : !! INPUTS
449 : : !! natom=number of atoms
450 : : !! nspden = number of spin densities (1 2 or 4)
451 : : !! option= if 0, the sum of the atomic spherical functions is returned in nv_constr_dft_r; if non-zero, they are added to nv_constr_dft_r
452 : : !! rprimd=lattice vectors (dimensionful)
453 : : !! mpi_enreg=mpi structure with communicator info
454 : : !! nfft=number of points in standard fft grid
455 : : !! ngfft=FFT grid dimensions
456 : : !! ntypat=number of types of atoms
457 : : !! ratsph(ntypat)=radii for muffin tin spheres of each atom
458 : : !! typat(natom)=types of atoms
459 : : !! xred(3,natom)=reduced atomic positions
460 : : !!
461 : : !! SIDE EFFECTS
462 : : !! nv_constr_dft_r=the constrained potential or density in real space
463 : : !!
464 : : !! SOURCE
465 : :
466 : 1403 : subroutine add_atomic_fcts(natom,nspden,rprimd,mpi_enreg,nfft,ngfft,ntypat,option,ratsph, &
467 : 1403 : ratsm, typat,coeffs_constr_dft,nv_constr_dft_r,xred)
468 : :
469 : : !Arguments ------------------------------------
470 : : !scalars
471 : : integer,intent(in) :: natom,nfft,nspden,ntypat,option
472 : : type(MPI_type),intent(in) :: mpi_enreg
473 : : !arrays
474 : : integer,intent(in) :: typat(natom)
475 : : integer,intent(in) :: ngfft(18)
476 : : real(dp),intent(in) :: coeffs_constr_dft(nspden,natom)
477 : : real(dp),intent(inout) :: nv_constr_dft_r(nfft,nspden)
478 : : real(dp),intent(in) :: ratsph(ntypat)
479 : : real(dp),intent(in) :: ratsm ! Ben change
480 : : real(dp),intent(in) :: rprimd(3,3)
481 : : real(dp),intent(in) :: xred(3,natom)
482 : :
483 : : !Local variables-------------------------------
484 : : !scalars
485 : : integer,parameter :: ishift=5
486 : : integer :: iatom, ierr
487 : : integer :: n1a, n1b, n3a, n3b, n2a, n2b
488 : : integer :: n1, n2, n3
489 : : integer :: ifft_local
490 : : integer :: i1,i2,i3,ix,iy,iz,izloc
491 : : real(dp) :: dfsm,dify,difz,fsm,r2atsph,rr1,rr2,rr3,ratsm2,rx23,ry23,rz23 !Ben change: remove ratsm
492 : : real(dp) :: r2,r2_11,r2_123,r2_23
493 : : real(dp) :: ucvol
494 : : real(dp),parameter :: delta=0.99_dp
495 : : !arrays
496 : 1403 : real(dp), allocatable :: difx(:)
497 : : real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3)
498 : : real(dp) :: tsec(2)
499 : 1403 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
500 : 1403 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
501 : : ! ***********************************************************************************************
502 : :
503 : : !We need the metric because it is needed to compute the "box" around each atom
504 : 1403 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
505 : :
506 : 1403 : n1 = ngfft(1)
507 : 1403 : n2 = ngfft(2)
508 : 1403 : n3 = ngfft(3)
509 : :
510 : : !Ben change: comment out ratsm line
511 : : !ratsm = 0.05_dp ! default value for the smearing region radius - may become input variable later
512 : :
513 [ + + ]: 1403 : if(option==0)then
514 [ + + + + ]: 1425865 : nv_constr_dft_r = zero
515 : : endif
516 : :
517 : : !Get the distrib associated with this fft_grid
518 : 1403 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
519 : :
520 : : !Loop over atoms
521 : : !-------------------------------------------
522 [ + + ]: 4634 : do iatom=1,natom
523 : :
524 [ + + + + ]: 14417 : if(sum(coeffs_constr_dft(1:nspden,iatom)**2)<tol12)then
525 : : cycle
526 : : endif
527 : :
528 : : ! Define a "box" around the atom
529 : 3005 : r2atsph=1.0000001_dp*ratsph(typat(iatom))**2
530 : 3005 : rr1=sqrt(r2atsph*gmet(1,1))
531 : 3005 : rr2=sqrt(r2atsph*gmet(2,2))
532 : 3005 : rr3=sqrt(r2atsph*gmet(3,3))
533 : :
534 : 3005 : n1a=int((xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1
535 : 3005 : n1b=int((xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1
536 : 3005 : n2a=int((xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2
537 : 3005 : n2b=int((xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2
538 : 3005 : n3a=int((xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3
539 : 3005 : n3b=int((xred(3,iatom)+rr3+ishift)*n3 )-ishift*n3
540 : :
541 : 3005 : ratsm2 = -(ratsm**2 - 2*ratsph(typat(iatom))*ratsm)
542 : :
543 [ - + + - : 9015 : ABI_MALLOC(difx,(n1a:n1b))
+ - + - -
+ - + ]
544 [ + + ]: 30074 : do i1=n1a,n1b
545 : 30074 : difx(i1)=dble(i1)/dble(n1)-xred(1,iatom)
546 : : enddo ! i1
547 : :
548 [ + + ]: 28936 : do i3=n3a,n3b
549 : 25931 : iz=mod(i3+ishift*n3,n3)
550 [ + - ]: 28936 : if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
551 : 25931 : izloc = ffti3_local(iz+1) - 1
552 : 25931 : difz=dble(i3)/dble(n3)-xred(3,iatom)
553 [ + + ]: 372416 : do i2=n2a,n2b
554 : 346485 : iy=mod(i2+ishift*n2,n2)
555 : 346485 : dify=dble(i2)/dble(n2)-xred(2,iatom)
556 : 346485 : rx23=dify*rprimd(1,2)+difz*rprimd(1,3)
557 : 346485 : ry23=dify*rprimd(2,2)+difz*rprimd(2,3)
558 : 346485 : rz23=dify*rprimd(3,2)+difz*rprimd(3,3)
559 : 346485 : r2_23=rx23**2+ry23**2+rz23**2
560 : 346485 : r2_11=rprimd(1,1)**2+rprimd(2,1)**2+rprimd(3,1)**2
561 : 346485 : r2_123=2*(rprimd(1,1)*rx23+rprimd(2,1)*ry23+rprimd(3,1)*rz23)
562 [ + + ]: 6459021 : do i1=n1a,n1b
563 : 6086605 : r2=(difx(i1)*r2_11+r2_123)*difx(i1)+r2_23
564 [ + + ]: 6086605 : if (r2 > r2atsph) cycle
565 : 3237143 : call radsmear(dfsm,fsm,r2,r2atsph,ratsm2)
566 : 3237143 : ix=mod(i1+ishift*n1,n1)
567 : : ! Identify the fft indexes of the rectangular grid around the atom
568 : 3237143 : ifft_local=1+ix+n1*(iy+n2*izloc)
569 [ + + ]: 10912967 : nv_constr_dft_r(ifft_local,1:nspden)=nv_constr_dft_r(ifft_local,1:nspden) + fsm*coeffs_constr_dft(1:nspden,iatom)
570 : :
571 : : end do ! i1
572 : : end do ! i2
573 : : end if ! if this is my fft slice
574 : : end do ! i3
575 [ - + ]: 4634 : ABI_FREE(difx)
576 : :
577 : : ! end loop over atoms
578 : : end do
579 : :
580 : : !MPI parallelization
581 : : !TODO: test if xmpi_sum does the correct stuff for a slice of nv_constr_dft_r
582 [ - + ]: 1403 : if(mpi_enreg%nproc_fft>1)then
583 : 0 : call timab(48,1,tsec)
584 : 0 : call xmpi_sum(nv_constr_dft_r,mpi_enreg%comm_fft,ierr)
585 : 0 : call timab(48,2,tsec)
586 : : end if
587 : :
588 : : ! write (201,*) '# potential 1'
589 : : ! write (201,*) nv_constr_dft_r(:,1)
590 : :
591 : : ! write (202,*) '# potential 2'
592 : : ! write (202,*) nv_constr_dft_r(:,2)
593 : :
594 : : ! if (nspden > 2) then
595 : : ! write (203,*) '# potential 3'
596 : : ! write (203,*) nv_constr_dft_r(:,3)
597 : :
598 : : ! write (204,*) '# potential 4'
599 : : ! write (204,*) nv_constr_dft_r(:,4)
600 : : ! end if
601 : :
602 : 1403 : end subroutine add_atomic_fcts
603 : : !!***
604 : :
605 : : !!****f* m_dens/constrained_dft_ini
606 : : !! NAME
607 : : !! constrained_dft_ini
608 : : !!
609 : : !! FUNCTION
610 : : !! Initialize the constrained_dft datastructure.
611 : : !! Mostly copying already available (dtset) information, but also computing intgf2
612 : :
613 : : !!
614 : : !! INPUTS
615 : : !! chrgat(natom) = target charge for each atom. Not always used, it depends on the value of constraint_kind
616 : : !! constraint_kind(ntypat)=for each type of atom, 0=no constraint,
617 : : !! 1=fix only the magnetization direction, following spinat direction,
618 : : !! 2=fix the magnetization vector to be the spinat one,
619 : : !! 3=fix the magnetization amplitude to be the spinat one, but does not fix its direction
620 : : !! other future values will constrain the local atomic charge and possibly mix constraints if needed.
621 : : !! magconon=type of penalty function (so, not constrained DFT).
622 : : !! magcon_lambda=strength of the atomic spherical constraint
623 : : !! mpi_enreg=mpi structure with communicator info
624 : : !! natom=number of atoms
625 : : !! nfft=number of points in standard fft grid
626 : : !! ngfft=FFT grid dimensions
627 : : !! nspden = number of spin densities (1 2 or 4)
628 : : !! ntypat=number of types of atoms
629 : : !! ratsm=smearing width for ratsph
630 : : !! ratsph(ntypat)=radii for muffin tin spheres of each atom
631 : : !! rprimd=lattice vectors (dimensioned)
632 : : !! spinat(3,natom)=magnetic moments vectors, possible targets according to the value of constraint_kind
633 : : !! typat(natom)=types of atoms
634 : : !! xred(3,natom)=reduced atomic positions
635 : : !! ziontypat(ntypat)=ionic charge, per type of atom
636 : : !!
637 : : !! OUTPUT
638 : : !! constrained_dft=datastructure that contain the needed information to enforce the density and magnetization constraints
639 : : !! Most of the data are simply copied from dtset, but also constrained_dft%intgf2(natom,natom) is computed from the available data.
640 : : !!
641 : : !! SOURCE
642 : :
643 [ - + - + : 64 : subroutine constrained_dft_ini(chrgat,constrained_dft,constraint_kind,&
- + - + -
+ - + -
+ ]
644 : : & magconon,magcon_lambda,mpi_enreg,natom,nfftf,ngfftf,nspden,ntypat,&
645 : 64 : & ratsm,ratsph,rprimd,spinat,typat,xred,ziontypat,znucl,qgbt,use_gbt,spinaxis)
646 : :
647 : : !Arguments ------------------------------------
648 : : !scalars
649 : : integer,intent(in) :: magconon,natom,nfftf,nspden,ntypat,use_gbt
650 : : real(dp),intent(in) :: magcon_lambda,ratsm
651 : : type(MPI_type),intent(in) :: mpi_enreg
652 : : type(constrained_dft_t),intent(out):: constrained_dft
653 : : !arrays
654 : : integer,intent(in) :: constraint_kind(ntypat)
655 : : integer,intent(in) :: ngfftf(18)
656 : : integer,intent(in) :: typat(natom)
657 : : real(dp),intent(in) :: chrgat(natom)
658 : : real(dp),intent(in) :: ratsph(ntypat)
659 : : real(dp),intent(in) :: rprimd(3,3)
660 : : real(dp),intent(in) :: spinat(3,natom)
661 : : real(dp),intent(in) :: xred(3,natom)
662 : : real(dp),intent(in) :: ziontypat(ntypat)
663 : : real(dp),intent(in) :: znucl(ntypat)
664 : : real(dp),intent(in) :: qgbt(3)
665 : : real(dp),intent(in) :: spinaxis(3)
666 : :
667 : : !Local variables-------------------------------
668 : : !scalars
669 : : integer,parameter :: cplex1=1
670 : : real(dp) :: ucvol
671 : : !arrays
672 : 64 : real(dp), allocatable :: intgf2(:,:) ! natom,natom
673 [ - + ]: 64 : real(dp), allocatable :: rhor_dum(:,:) ! nfftf,nspden
674 : : real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3)
675 : : ! ***********************************************************************************************
676 : :
677 [ - + + - : 256 : ABI_MALLOC(intgf2,(natom,natom))
- + + - +
- + - - +
- + ]
678 : :
679 : : !We need the metric because it is needed in calcdenmagsph.F90
680 : 64 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
681 : :
682 [ + + + + : 75 : if(any(constraint_kind(:)/=0))then
+ + ]
683 : : !We need to precompute intgf2
684 [ - + + - : 216 : ABI_MALLOC(rhor_dum,(nfftf,nspden))
- + + - +
- + - - +
- + ]
685 [ + + + + ]: 2489335 : rhor_dum(:,:)=zero
686 : : call calcdenmagsph(mpi_enreg,natom,nfftf,ngfftf,nspden,ntypat,&
687 : 54 : & ratsm,ratsph,rhor_dum,rprimd,typat,xred,0,cplex1,qgbt,use_gbt,intgf2=intgf2)
688 [ - + ]: 54 : ABI_FREE(rhor_dum)
689 : : else
690 [ + + + + ]: 30 : intgf2=zero
691 : : endif
692 : :
693 : 64 : constrained_dft%magconon =magconon
694 : 64 : constrained_dft%magcon_lambda =magcon_lambda
695 : 64 : constrained_dft%natom =natom
696 : 64 : constrained_dft%nfftf =nfftf
697 [ + + ]: 1216 : constrained_dft%ngfftf =ngfftf
698 : 64 : constrained_dft%nspden =nspden
699 : 64 : constrained_dft%ntypat =ntypat
700 : 64 : constrained_dft%ratsm =ratsm
701 [ + + + + ]: 832 : constrained_dft%rprimd =rprimd
702 : :
703 [ - + + - : 192 : ABI_MALLOC(constrained_dft%chrgat,(natom))
+ - + - -
+ - + -
+ ]
704 [ - + + - : 192 : ABI_MALLOC(constrained_dft%constraint_kind,(ntypat))
+ - + - -
+ - + -
+ ]
705 [ - + + - : 192 : ABI_MALLOC(constrained_dft%intgf2,(natom,natom))
- + + - -
+ - + -
+ ]
706 [ - + + - : 192 : ABI_MALLOC(constrained_dft%ratsph,(ntypat))
+ - + - -
+ - + -
+ ]
707 [ - + + - : 192 : ABI_MALLOC(constrained_dft%spinat,(3,natom))
+ - + - -
+ - + -
+ ]
708 [ - + + - : 192 : ABI_MALLOC(constrained_dft%typat,(natom))
+ - + - -
+ - + -
+ ]
709 [ - + + - : 128 : ABI_MALLOC(constrained_dft%ziontypat,(ntypat))
- + - + -
+ ]
710 [ - + + - : 128 : ABI_MALLOC(constrained_dft%znucl,(ntypat))
- + - + -
+ ]
711 : :
712 [ - + + + ]: 264 : constrained_dft%chrgat =chrgat
713 [ - + + + ]: 209 : constrained_dft%constraint_kind =constraint_kind
714 [ - + + + : 592 : constrained_dft%intgf2 =intgf2
+ + ]
715 [ - + + + ]: 209 : constrained_dft%ratsph =ratsph
716 [ - + + + : 672 : constrained_dft%spinat =spinat
+ + ]
717 [ - + + + ]: 264 : constrained_dft%typat =typat
718 [ - + + + ]: 209 : constrained_dft%ziontypat =ziontypat
719 [ - + + + ]: 209 : constrained_dft%znucl =znucl
720 [ + + ]: 256 : constrained_dft%spinaxis =spinaxis
721 : :
722 [ - + ]: 64 : ABI_FREE(intgf2)
723 : :
724 : 64 : end subroutine constrained_dft_ini
725 : : !!***
726 : :
727 : :
728 : : !!****f* m_dens/constrained_dft_free
729 : : !! NAME
730 : : !! constrained_dft_free
731 : : !!
732 : : !! FUNCTION
733 : : !! Free the constrained_dft datastructure.
734 : : !!
735 : : !! SOURCE
736 : :
737 : 6962 : subroutine constrained_dft_free(constrained_dft)
738 : :
739 : : !Arguments ------------------------------------
740 : : type(constrained_dft_t),intent(inout):: constrained_dft
741 : : ! ***********************************************************************************************
742 : :
743 [ + + ]: 6962 : ABI_SFREE(constrained_dft%chrgat)
744 [ + + ]: 6962 : ABI_SFREE(constrained_dft%constraint_kind)
745 [ + + ]: 6962 : ABI_SFREE(constrained_dft%intgf2)
746 [ + + ]: 6962 : ABI_SFREE(constrained_dft%ratsph)
747 [ + + ]: 6962 : ABI_SFREE(constrained_dft%spinat)
748 [ + + ]: 6962 : ABI_SFREE(constrained_dft%typat)
749 [ + + ]: 6962 : ABI_SFREE(constrained_dft%ziontypat)
750 [ + + ]: 6962 : ABI_SFREE(constrained_dft%znucl)
751 : :
752 : 6962 : end subroutine constrained_dft_free
753 : : !!***
754 : :
755 : :
756 : : !!****f* m_dens/constrained_residual
757 : : !! NAME
758 : : !! constrained_residual
759 : : !!
760 : : !! FUNCTION
761 : : !! Recompute the residual to take into account the constraints, within constrained DFT.
762 : : !! The kind of constraint is given by constraint_kind, and the target values are given by spinat, for the local atomic magnetization,
763 : : !! and chrgat minus the ionic charge, for the local atomic charge.
764 : :
765 : : !!
766 : : !! INPUTS
767 : : !! c_dft <type(constrained_dft_t)>=datastructure for the information related to constrained DFT
768 : : !! ! chrgat(natom) = target charge for each atom. Not always used, it depends on the value of constraint_kind
769 : : !! ! constraint_kind(ntypat)=for each type of atom, 0=no constraint,
770 : : !! ! 1=fix only the magnetization direction, following spinat direction,
771 : : !! ! 2=fix the magnetization vector to be the spinat one,
772 : : !! ! 3=fix the magnetization amplitude to be the spinat one, but does not fix its direction
773 : : !! ! other future values will constrain the local atomic charge and possibly mix constraints if needed.
774 : : !! ! intgf2(natom,natom)=(precomputed) overlap of the spherical integration functions for each atom in a sphere of radius ratsph.
775 : : !! ! magcon_lambda=strength of the atomic spherical constraint
776 : : !! ! natom=number of atoms
777 : : !! ! nfftf=number of points in fine fft grid
778 : : !! ! ngfftf=FFT grid dimensions
779 : : !! ! nspden = number of spin densities (1 2 or 4)
780 : : !! ! ntypat=number of types of atoms
781 : : !! ! ratsm=smearing width for ratsph
782 : : !! ! ratsph(ntypat)=radii for muffin tin spheres of each atom
783 : : !! ! rprimd=lattice vectors (dimensioned)
784 : : !! ! spinat(3,natom)=magnetic moments vectors, possible targets according to the value of constraint_kind
785 : : !! ! typat(natom)=types of atoms
786 : : !! ! ziontypat(ntypat)= ionic charge, per type of atom
787 : : !! mpi_enreg=mpi structure with communicator info
788 : : !! rhor(nfft,nspden)=array for electron density in el./bohr**3. At output it will be constrained.
789 : : !! xred(3,natom)=reduced atomic positions
790 : : !!
791 : : !! OUTPUT
792 : : !! e_constrained_dft=correction to the total energy, to make it variational
793 : : !! grcondft(3,natom)=d(E_constrained_DFT)/d(xred) (hartree)
794 : : !! intgres(nspden,natom)=integrated residuals from constrained DFT. They are also Lagrange parameters, or gradients with respect to constraints.
795 : : !! strscondft(6)=stress due to constraints = -d(E_constrained_DFT)/d(strain) / ucvol (hartree/Bohr^3)
796 : : !!
797 : : !! SIDE EFFECTS
798 : : !! vresid(nfft,nspden)==array for potential residual in real space
799 : : !! At output it will be modified: projected onto the space orthogonal to the atomic spherical functions (if there is a related
800 : : !! constrained, and augmented by such atomic spherical functions multiplied by the difference between the actual
801 : : !! integrated charge or magnetization and the target ones.
802 : : !!
803 : : !! SOURCE
804 : :
805 [ - + - + ]: 1294 : subroutine constrained_residual(c_dft,e_constrained_dft,grcondft,intgres,mpi_enreg,rhor,strscondft,vresid,xred,qgbt,use_gbt)
806 : :
807 : : !Arguments ------------------------------------
808 : : !scalars
809 : : class(constrained_dft_t),intent(in) :: c_dft
810 : : integer, intent(in) :: use_gbt
811 : : real(dp),intent(out) :: e_constrained_dft
812 : : type(MPI_type),intent(in) :: mpi_enreg
813 : : !arrays
814 : : real(dp),intent(out) :: grcondft(:,:) ! 3,natom
815 : : real(dp),intent(out) :: intgres(:,:) ! nspden,natom
816 : : real(dp),intent(in) :: rhor(c_dft%nfftf,c_dft%nspden)
817 : : real(dp),intent(out) :: strscondft(6)
818 : : real(dp),intent(inout) :: vresid(c_dft%nfftf,c_dft%nspden)
819 : : real(dp),intent(in) :: xred(3,c_dft%natom),qgbt(3)
820 : :
821 : : !Local variables-------------------------------
822 : : !scalars
823 : : integer :: conkind,iatom,ii,jatom,info,natom,nfftf,nspden,ntypat,option
824 : : integer,parameter :: cplex1=1
825 : : real(dp) :: intgd,intgden_norm,intgden_proj,intgres_proj,norm,scprod
826 : : !arrays
827 : 2588 : integer :: ipiv(c_dft%natom)
828 : : real(dp) :: corr_denmag(4),gr_intgd(3),strs_intgd(6)
829 : 1294 : real(dp), allocatable :: coeffs_constr_dft(:,:) ! nspden,natom
830 : 1294 : real(dp), allocatable :: gr_intgden(:,:,:) ! 3,nspden,natom
831 : 1294 : real(dp), allocatable :: intgden(:,:) ! nspden,natom
832 : 1294 : real(dp), allocatable :: intgden_delta(:,:) ! nspden,natom
833 : 1294 : real(dp), allocatable :: intgres_tmp(:,:) ! nspden,natom
834 : 1294 : real(dp), allocatable :: intgr(:,:) ! natom,nspden
835 : 1294 : real(dp), allocatable :: strs_intgden(:,:,:) ! 6,nspden,natom
836 : 2588 : real(dp) :: intgf2(c_dft%natom,c_dft%natom),rhomag(2,c_dft%nspden),work(2*c_dft%natom)
837 : : real(dp) :: intgden_normed(3)
838 : : real(dp) :: spinat_normed(3)
839 : : ! ***********************************************************************************************
840 : :
841 : : !DEBUG
842 : : !write(std_out,*) ' constrained_residual : enter '
843 : : !ENDDEBUG
844 : :
845 : 1294 : natom=c_dft%natom
846 : 1294 : nfftf=c_dft%nfftf
847 : 1294 : nspden=c_dft%nspden
848 : 1294 : ntypat=c_dft%ntypat
849 : :
850 : : !We need the integrated magnetic moments
851 [ - + + - : 5176 : ABI_MALLOC(intgden,(nspden,natom))
- + + - +
- + - - +
- + ]
852 [ - + + - : 5176 : ABI_MALLOC(gr_intgden,(3,nspden,natom))
- + + - +
- + - - +
- + ]
853 [ - + + - : 5176 : ABI_MALLOC(strs_intgden,(6,nspden,natom))
- + + - +
- + - - +
- + ]
854 : :
855 : : call calcdenmagsph(mpi_enreg,natom,nfftf,c_dft%ngfftf,nspden,ntypat,c_dft%ratsm,c_dft%ratsph,rhor,c_dft%rprimd,c_dft%typat,&
856 : 1294 : xred,1,cplex1,qgbt,use_gbt,intgden=intgden,gr_intgden=gr_intgden,rhomag=rhomag,strs_intgden=strs_intgden)
857 [ + + ]: 2588 : call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,[std_out],1,qgbt,c_dft%ratsm,c_dft%ratsph,rhomag,c_dft%typat,c_dft%znucl,c_dft%spinaxis)
858 : :
859 : : !DEBUG
860 : : !write(std_out,*) ' intgden(1:nspden,1:natom)=',intgden(1:nspden,1:natom)
861 : : !ENDDEBUG
862 : :
863 : : !We need the integrated residuals
864 [ - + + - : 3882 : ABI_MALLOC(intgres_tmp,(nspden,natom))
- + + - -
+ - + ]
865 [ + + + + ]: 15254 : intgres_tmp(:,:)=zero
866 : : call calcdenmagsph(mpi_enreg,natom,nfftf,c_dft%ngfftf,nspden,ntypat,&
867 : 1294 : & c_dft%ratsm,c_dft%ratsph,vresid,c_dft%rprimd,c_dft%typat,xred,11,cplex1,qgbt,use_gbt,intgden=intgres_tmp,rhomag=rhomag)
868 : :
869 : : !DEBUG
870 : : !write(std_out,*) ' intgres_tmp(1:nspden,1:natom)=',intgres_tmp(1:nspden,1:natom)
871 : : !ENDDEBUG
872 : :
873 : : !Make the proper combination of intgres_tmp, to single out the scalar potential residual and the magnetic field potential residuals for x,y,z.
874 [ + + + + ]: 15254 : intgres(:,:)=zero
875 [ + + ]: 4416 : do iatom=1,natom
876 [ + + ]: 3122 : if(nspden==1)then
877 : 210 : intgres(1,iatom)=intgres_tmp(1,iatom)
878 [ + + ]: 2912 : else if(nspden==2)then
879 : 510 : intgres(1,iatom)=half*(intgres_tmp(1,iatom)+intgres_tmp(2,iatom))
880 : 510 : intgres(2,iatom)=half*(intgres_tmp(1,iatom)-intgres_tmp(2,iatom))
881 [ + - ]: 2402 : else if(nspden==4)then
882 : : !Change the potential residual to the density+magnetization convention
883 : 2402 : intgres(1,iatom)=half*(intgres_tmp(1,iatom)+intgres_tmp(2,iatom))
884 : 2402 : intgres(2,iatom)= intgres_tmp(3,iatom)
885 : 2402 : intgres(3,iatom)=-intgres_tmp(4,iatom)
886 : 2402 : intgres(4,iatom)=half*(intgres_tmp(1,iatom)-intgres_tmp(2,iatom))
887 : : endif
888 : 3122 : conkind=c_dft%constraint_kind(c_dft%typat(iatom))
889 [ + + ]: 3122 : if(conkind <10)intgres(1,iatom)=zero
890 [ + + + + : 4641 : if( mod(conkind,10)==0 .and. nspden>1)intgres(2:nspden,iatom)=zero
+ + ]
891 : : enddo
892 : : !Print the potential residuals
893 [ - + - - : 2588 : call prtdenmagsph(cplex1,intgres,natom,nspden,ntypat,[std_out],11,qgbt,c_dft%ratsm,c_dft%ratsph,rhomag,c_dft%typat,c_dft%znucl,c_dft%spinaxis)
- - - - +
+ - + ]
894 [ - + ]: 1294 : ABI_FREE(intgres_tmp)
895 : :
896 : : !Also exchanges the spin and atom indices to prepare the solution of the linear system of equation
897 [ - + + - : 3882 : ABI_MALLOC(intgr,(natom,nspden))
- + + - -
+ - + ]
898 [ + + ]: 4416 : do iatom=1,natom
899 [ + + ]: 15254 : intgr(iatom,1:nspden)=intgres(1:nspden,iatom)
900 : : enddo
901 : :
902 : : !In case there is an overlap between spheres, must solve the linear system of equations
903 : : !and take into account non-diagonal elements only for the set of atoms for which there is a constraint.
904 : : !This can be different for the
905 : : !charge residual or for the spin residual (but for the spin constraints, the set of atoms is inclusive of all constraints)
906 [ + + ]: 3882 : do ii=1,2 ! Charge, then spin
907 [ + + + + ]: 25308 : intgf2(:,:)=zero
908 [ + + ]: 8832 : do iatom=1,natom
909 : 8832 : intgf2(iatom,iatom)=c_dft%intgf2(iatom,iatom)
910 : : enddo
911 [ + + ]: 6244 : do iatom=1,natom-1
912 [ + + ]: 15016 : do jatom=iatom,natom
913 [ + + ]: 12428 : if(c_dft%intgf2(iatom,jatom)>tol8)then
914 : : !In the charge case, must have both atoms with constraints bigger than 10
915 [ + + ]: 5880 : if(ii==1)then
916 [ + + + + ]: 2940 : if(c_dft%constraint_kind(c_dft%typat(iatom))>=10 .and. &
917 : : & c_dft%constraint_kind(c_dft%typat(jatom))>=10 )then
918 : 609 : intgf2(iatom,jatom)=c_dft%intgf2(iatom,jatom)
919 : 609 : intgf2(jatom,iatom)=c_dft%intgf2(iatom,jatom)
920 : : endif
921 : : endif
922 : : !In the spin case, must have both atoms with constraints not ending with 0
923 [ + + ]: 5880 : if(ii==2)then
924 [ + + + + ]: 2940 : if(mod(c_dft%constraint_kind(c_dft%typat(iatom)),10)/=0 .and. &
925 : : & mod(c_dft%constraint_kind(c_dft%typat(jatom)),10)/=0 )then
926 : 2728 : intgf2(iatom,jatom)=c_dft%intgf2(iatom,jatom)
927 : 2728 : intgf2(jatom,iatom)=c_dft%intgf2(iatom,jatom)
928 : : endif
929 : : endif
930 : : endif
931 : : enddo
932 : : enddo
933 : :
934 : : !Solve the linear system of equation, for the different spins
935 : 2588 : call dsytrf('U',natom,intgf2,natom,ipiv,work,2*natom,info)
936 : : ! call dsytri('U',natom,intgf2,natom,ipiv,work,info)
937 [ + + ]: 2588 : if(ii==1)then
938 : 1294 : call dsytrs('U',natom,1,intgf2,natom,ipiv,intgr,natom,info)
939 [ + + ]: 1294 : else if(ii==2 .and. nspden>1)then
940 [ - + - - : 1189 : call dsytrs('U',natom,nspden-1,intgf2,natom,ipiv,intgr(1:natom,2:nspden),natom,info)
- - - - -
+ - - -
- ]
941 : : endif
942 : :
943 : : !Store the new residuals
944 [ + + ]: 10126 : do iatom=1,natom
945 [ + + ]: 6244 : if(ii==1)intgres(1,iatom)=intgr(iatom,1)
946 [ + + + + ]: 16548 : if(ii==2)intgres(2:nspden,iatom)=intgr(iatom,2:nspden)
947 : : enddo
948 : :
949 : : enddo
950 : :
951 : : !DEBUG
952 : : !write(std_out,*) ' after multiplication by ftt-1 , so, torque :'
953 : : !write(std_out,*) ' intgres(1:nspden,1:natom)=',intgres(1:nspden,1:natom)
954 : : !ENDDEBUG
955 : :
956 : : !Compute the delta of the integrated dens with respect to the target
957 : : !Compute the energy correction, to make the energy functional variational
958 : : !Also projects the residual in case constraint_kind 2
959 : 1294 : e_constrained_dft=zero
960 [ + + + + ]: 13782 : grcondft=zero
961 : 1294 : strscondft=zero
962 [ - + + - : 5176 : ABI_MALLOC(intgden_delta,(nspden,natom))
- + + - +
- + - - +
- + ]
963 [ + + + + ]: 15254 : intgden_delta(:,:)=zero
964 [ + + ]: 4416 : do iatom=1,natom
965 : :
966 : : ! The integrated density must be in the total density+magnetization representation
967 [ + + ]: 3122 : if(nspden==2)then
968 : 510 : intgd =intgden(1,iatom)+intgden(2,iatom)
969 : 510 : intgden(2,iatom)=intgden(1,iatom)-intgden(2,iatom)
970 : 510 : intgden(1,iatom)=intgd
971 [ + + ]: 2040 : do ii=1,3
972 : 1530 : gr_intgd(ii) =gr_intgden(ii,1,iatom)+gr_intgden(ii,2,iatom)
973 : 1530 : gr_intgden(ii,2,iatom)=gr_intgden(ii,1,iatom)-gr_intgden(ii,2,iatom)
974 : 2040 : gr_intgden(ii,1,iatom)=gr_intgd(ii)
975 : : enddo
976 [ + + ]: 3570 : do ii=1,6
977 : 3060 : strs_intgd(ii) =strs_intgden(ii,1,iatom)+strs_intgden(ii,2,iatom)
978 : 3060 : strs_intgden(ii,2,iatom)=strs_intgden(ii,1,iatom)-strs_intgden(ii,2,iatom)
979 : 3570 : strs_intgden(ii,1,iatom)=strs_intgd(ii)
980 : : enddo
981 : : endif
982 : :
983 : : !Comparison with the target value, and computation of the correction in terms of density and magnetization coefficients.
984 : 3122 : conkind=c_dft%constraint_kind(c_dft%typat(iatom))
985 : :
986 [ + + ]: 3122 : if(conkind >=10)then
987 : : !The electronic constraint is such that the ziontypat charge minus (the electronic charge is negative) the atomic electronic density
988 : : !intgden gives the target charge chrgat.
989 : 825 : intgden_delta(1,iatom)=intgden(1,iatom)+c_dft%chrgat(iatom)-c_dft%ziontypat(c_dft%typat(iatom))
990 : : ! Uses the usual electronic charge definition, instead of the total nucleus-electronic charge
991 : : ! intgden_delta(1,iatom)=intgden(1,iatom)-c_dft%chrgat(iatom)
992 : : endif
993 : :
994 [ + + + - ]: 3122 : if( mod(conkind,10)==1 .and. nspden>1)then
995 : : !Fix the different components of the magnetization vector
996 [ + + ]: 2438 : if(nspden==2)intgden_delta(2,iatom)=intgden(2,iatom)-c_dft%spinat(3,iatom)
997 [ + + + + ]: 8828 : if(nspden==4)intgden_delta(2:4,iatom)=intgden(2:4,iatom)-c_dft%spinat(1:3,iatom)
998 [ + + + - ]: 684 : else if( ( mod(conkind,10)>=2 .and. mod(conkind,10)<=4) .and. nspden>1)then
999 [ + + ]: 1372 : norm = sqrt(sum(c_dft%spinat(:,iatom)**2))
1000 [ + - ]: 343 : if (norm > tol10) then
1001 [ + + ]: 343 : if( mod(conkind,10)==2 )then
1002 : : !Fix the axis of the magnetization vector
1003 [ + + ]: 149 : if(nspden==4)then
1004 [ + + ]: 348 : spinat_normed(:) = c_dft%spinat(:,iatom) / norm
1005 : : !Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector
1006 : : !This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector
1007 : : intgden_proj=spinat_normed(1)*intgden(2,iatom)+ &
1008 : : & spinat_normed(2)*intgden(3,iatom)+ &
1009 : 87 : & spinat_normed(3)*intgden(4,iatom)
1010 [ + + ]: 348 : intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom)-spinat_normed(1:3)*intgden_proj
1011 : : !Also projects the residual, so that the usual optimization is done is the largest possible space
1012 : : intgres_proj=spinat_normed(1)*intgres(2,iatom)+ &
1013 : : & spinat_normed(2)*intgres(3,iatom)+ &
1014 : 87 : & spinat_normed(3)*intgres(4,iatom)
1015 [ + + ]: 348 : intgres(2:nspden,iatom)=intgres(2:nspden,iatom)-spinat_normed(1:3)*intgres_proj
1016 [ + - ]: 62 : else if(nspden==2)then
1017 : : !The direction must be correct, collinear, so no change.
1018 : 62 : intgden_delta(2,iatom)=zero
1019 : : endif
1020 [ + + ]: 194 : else if( mod(conkind,10)==3 )then
1021 : : !Fix the amplitude of the magnetization vector
1022 [ + + ]: 320 : intgden_norm = sqrt(sum(intgden(2:nspden,iatom)**2))
1023 [ + + ]: 320 : intgden_delta(2:nspden,iatom)=(one-norm/intgden_norm)*intgden(2:nspden,iatom)
1024 : : else if( mod(conkind,10)==4 )then
1025 : : !Fix the direction of the magnetization vector
1026 [ + - ]: 86 : if(nspden==4)then
1027 [ + + ]: 344 : spinat_normed(:) = c_dft%spinat(:,iatom) / norm
1028 [ + + ]: 344 : intgden_norm = sqrt(sum(intgden(2:nspden,iatom)**2))
1029 [ + + ]: 344 : intgden_normed(:) = intgden(2:nspden,iatom)/intgden_norm
1030 : : !Calculate the difference vector between the actual magnetization vectoro, times the scalar product between the
1031 : : !directions of present magnetization and target one, and the target magnetization direction renormalized by the magnetization length.
1032 : : !See notes 12 October 2021
1033 : : !DEBUG First possibility
1034 : : ! intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom)-spinat_normed(1:3)*intgden_norm
1035 : : !ENDDEBUG
1036 : : !DEBUG Second possibility
1037 [ + + ]: 344 : scprod=sum(spinat_normed(:)*intgden_normed(:))
1038 [ + + ]: 344 : intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom)*scprod-spinat_normed(1:3)*intgden_norm
1039 : : !ENDDEBUG
1040 : : !Also projects the residual. There might be some misalignement of the intgden_delta with the correct space of allowed variations
1041 : : !(not exactly perpendicular to spinat_normed or intgden_normed, or even frankly not at all perpendicular),
1042 : : !but when close to fulfilling the constraint, this difference becomes negligible.
1043 : : !The present choice, couple with the above definition of intgden_delta aligns both.
1044 : : !DEBUG First possibility
1045 : : intgres_proj=intgden_normed(1)*intgres(2,iatom)+ &
1046 : : & intgden_normed(2)*intgres(3,iatom)+ &
1047 : 86 : & intgden_normed(3)*intgres(4,iatom)
1048 [ + + ]: 344 : intgres(2:nspden,iatom)=intgres(2:nspden,iatom)-intgden_normed(1:3)*intgres_proj
1049 : : !ENDDEBUG
1050 : : !DEBUG Second possibility
1051 : : ! intgres_proj=spinat_normed(1)*intgres(2,iatom)+ &
1052 : : !& spinat_normed(2)*intgres(3,iatom)+ &
1053 : : !& spinat_normed(3)*intgres(4,iatom)
1054 : : ! intgres(2:nspden,iatom)=intgres(2:nspden,iatom)-spinat_normed(1:3)*intgres_proj
1055 : : !ENDDEBUG
1056 [ # # ]: 0 : else if(nspden==2)then
1057 : : !This case is ill-defined ... What follows is rather arbitrary.
1058 : : !When the actual magnetization and the target magnetization have same sign, the intgden_delta is zero.
1059 : : ! Otherwise, intgres is set to the difference between the present intgden and the target
1060 : 0 : intgden_delta(2,iatom)=zero
1061 [ # # ]: 0 : if(c_dft%spinat(2,iatom)*intgden(2,iatom)<-tol10)intgden_delta(2,iatom)=c_dft%spinat(2,iatom)-intgden(2,iatom)
1062 : : endif
1063 : : endif
1064 : : else
1065 : : !In this case (norm of constraint vanishes), we set the atomic magnetization to zero.
1066 [ # # ]: 0 : intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom)
1067 : : endif
1068 : : end if
1069 : : ! Lagrange energy contribution. Note that intgres is the derivative with respect to the chrgat and spinat constraints,
1070 : : ! because chrgat and spinat have directly been used in the definition of intgden_delta.
1071 : : ! WOOPS : chrgat comes with a POSITIVE sign in intgden_delta ?!?!
1072 [ + + ]: 13960 : e_constrained_dft=e_constrained_dft-sum(intgden_delta(:,iatom)*intgres(:,iatom))
1073 [ + + ]: 12488 : do ii=1,3
1074 [ + + ]: 45002 : grcondft(ii,iatom)=grcondft(ii,iatom)-sum(gr_intgden(ii,:,iatom)*intgres(:,iatom))
1075 : : enddo
1076 : : ! For the stress, this is the place where the summation over atoms is performed.
1077 [ + + ]: 23148 : do ii=1,6
1078 [ + + ]: 86882 : strscondft(ii)=strscondft(ii)-sum(strs_intgden(ii,:,iatom)*intgres(:,iatom))
1079 : : enddo
1080 : :
1081 : : !DEBUG
1082 : : ! write(6,*)' calcdenmagsph/constrained_residual, line 1058 : iatom=',iatom
1083 : : ! write(6,*)' e_constrained_dft,intgden_delta(:,iatom),intgres(:,iatom)=',e_constrained_dft,intgden_delta(:,iatom),intgres(:,iatom)
1084 : : ! write(6,*)' grcondft(1,iatom),gr_intgden(1,:,iatom),intgres(:,iatom)=',grcondft(1,iatom),gr_intgden(1,:,iatom),intgres(:,iatom)
1085 : : ! write(6,*)' strscondft(1),strs_intgden(1,:,iatom),intgres(:,iatom)=',strscondft(1),strs_intgden(1,:,iatom),intgres(:,iatom)
1086 : : ! call flush(6)
1087 : : ! stop
1088 : : !ENDDEBUG
1089 : :
1090 : : enddo
1091 : :
1092 [ - + ]: 1294 : ABI_FREE(gr_intgden)
1093 [ - + ]: 1294 : ABI_FREE(strs_intgden)
1094 : :
1095 [ - + + - : 3882 : ABI_MALLOC(coeffs_constr_dft,(nspden,natom))
- + + - -
+ - + ]
1096 [ + + + + ]: 15254 : coeffs_constr_dft=zero
1097 : :
1098 : : !With the delta of the integrated density and the atomic residual, compute the atomic correction to be applied to the potential
1099 : : !See the Eqs.(31), (33) and (34) of notes.
1100 [ + + ]: 4416 : do iatom=1,natom
1101 : :
1102 : : !Computation of the correction in terms of density and magnetization coefficients.
1103 : 3122 : conkind=c_dft%constraint_kind(c_dft%typat(iatom))
1104 : 3122 : corr_denmag(:)=zero
1105 : :
1106 [ + + ]: 3122 : if(conkind >=10)then
1107 : 825 : corr_denmag(1)=intgden_delta(1,iatom)*c_dft%magcon_lambda - intgres(1,iatom)
1108 : : endif
1109 : :
1110 [ + + + - ]: 3122 : if( mod(conkind,10)==1 .and. nspden>1)then
1111 : :
1112 : : !Fix the different components of the magnetization vector
1113 [ + + ]: 2438 : if(nspden==2)corr_denmag(2)=intgden_delta(2,iatom)*c_dft%magcon_lambda - intgres(2,iatom)
1114 [ + + + + ]: 8828 : if(nspden==4)corr_denmag(2:4)=intgden_delta(2:4,iatom)*c_dft%magcon_lambda - intgres(2:4,iatom)
1115 : :
1116 [ + + + - ]: 684 : else if( ( mod(conkind,10)>=2 .and. mod(conkind,10)<=4) .and. nspden>1)then
1117 : :
1118 [ + + ]: 1372 : norm = sqrt(sum(c_dft%spinat(:,iatom)**2))
1119 [ + - ]: 343 : if (norm > tol10) then
1120 : :
1121 [ + + ]: 343 : if( mod(conkind,10)==2 .or. mod(conkind,10)==4)then
1122 : : !Fix the axis/direction of the magnetization vector
1123 [ + + - + ]: 235 : if(nspden==4 .or. mod(conkind,10)==4)then
1124 [ + + ]: 692 : corr_denmag(2:nspden)=intgden_delta(2:nspden,iatom)*c_dft%magcon_lambda -intgres(2:nspden,iatom)
1125 : : else if(nspden==2 .and. mod(conkind,10)==2)then
1126 : : !The direction must be correct, collinear, so no change.
1127 : : corr_denmag(2)=zero
1128 : : endif
1129 : :
1130 : : else if( mod(conkind,10)==3 )then
1131 : :
1132 : : !Fix the amplitude of the magnetization vector
1133 : : !This is a special case, one does not work (at present) with the intgden_delta, while one should ...
1134 [ + + ]: 320 : intgden(2:nspden,iatom)=intgden(2:nspden,iatom) - intgres(2:nspden,iatom)/c_dft%magcon_lambda
1135 [ + + ]: 320 : intgden_norm = sqrt(sum(intgden(2:nspden,iatom)**2))
1136 [ + + ]: 320 : corr_denmag(2:nspden)=(one-norm/intgden_norm)*intgden(2:nspden,iatom)
1137 [ + + ]: 320 : corr_denmag(2:nspden)=corr_denmag(2:nspden)*c_dft%magcon_lambda
1138 : :
1139 : : endif
1140 : :
1141 : : else
1142 : : !In this case, we set the atomic magnetization to zero.
1143 [ # # ]: 0 : corr_denmag(2:nspden)=intgden_delta(2:nspden,iatom)*c_dft%magcon_lambda - intgres(2,iatom)
1144 : : endif
1145 : :
1146 : : end if
1147 : :
1148 : : !Convert from density/magnetization constraint residual to actual coefficient that will multiply the spherical function for the potential
1149 [ + + ]: 4416 : if(nspden==1)then
1150 : : !From charge to potential
1151 : 210 : coeffs_constr_dft(1,iatom)=corr_denmag(1)
1152 [ + - ]: 2912 : else if(nspden==2 .or. nspden==4)then
1153 : : !From charge and magnetization to potential
1154 : 2912 : coeffs_constr_dft(1,iatom)=corr_denmag(1)+corr_denmag(nspden)
1155 : 2912 : coeffs_constr_dft(2,iatom)=corr_denmag(1)-corr_denmag(nspden)
1156 [ + + ]: 2912 : if(nspden==4)then
1157 : 2402 : coeffs_constr_dft(3,iatom)= corr_denmag(2)
1158 : 2402 : coeffs_constr_dft(4,iatom)=-corr_denmag(3)
1159 : : endif
1160 : : endif
1161 : :
1162 : : enddo
1163 : :
1164 : : !Now compute the new residual, by adding the spherical functions
1165 : 1294 : option=1
1166 : : call add_atomic_fcts(natom,nspden,c_dft%rprimd,mpi_enreg,nfftf,c_dft%ngfftf,ntypat,option,&
1167 : 1294 : & c_dft%ratsph,c_dft%ratsm,c_dft%typat,coeffs_constr_dft,vresid,xred) ! Ben change: add ratsm
1168 : :
1169 [ - + ]: 1294 : ABI_FREE(coeffs_constr_dft)
1170 [ - + ]: 1294 : ABI_FREE(intgden)
1171 [ - + ]: 1294 : ABI_FREE(intgden_delta)
1172 [ - + ]: 1294 : ABI_FREE(intgr)
1173 : :
1174 : 1294 : end subroutine constrained_residual
1175 : : !!***
1176 : :
1177 : : !!****f* m_dens/mag_penalty
1178 : : !! NAME
1179 : : !! mag_penalty
1180 : : !!
1181 : : !! FUNCTION
1182 : : !! This routine is called to compute the potential corresponding to constrained magnetic moments using the penalty function algorithm.
1183 : : !!
1184 : : !! INPUTS
1185 : : !! c_dft <type(constrained_dft_t)>=datastructure for the information related to constrained DFT
1186 : : !! ! magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size
1187 : : !! ! magcon_lambda=strength of the atomic spherical constraint
1188 : : !! ! natom=number of atoms
1189 : : !! ! nfftf=number of points in fine fft grid
1190 : : !! ! ngfftf=FFT grid dimensions
1191 : : !! ! nspden = number of spin densities (1 2 or 4)
1192 : : !! ! ntypat=number of types of atoms
1193 : : !! ! ratsm=smearing width for ratsph
1194 : : !! ! ratsph(ntypat)=radii for muffin tin spheres of each atom
1195 : : !! ! rprimd=lattice vectors (dimensioned)
1196 : : !! ! spinat(3,natom)=magnetic moments vectors, possible targets according to the value of constraint_kind
1197 : : !! ! typat(natom)=types of atoms
1198 : : !! mpi_enreg=mpi structure with communicator info
1199 : : !! rhor=density in real space
1200 : : !! xred=reduced atomic positions
1201 : : !!
1202 : : !! OUTPUT
1203 : : !! nv_constr_dft_r=the constrained potential
1204 : : !!
1205 : : !! NOTES
1206 : : !! based on html notes for the VASP implementation at
1207 : : !! http://cms.mpi.univie.ac.at/vasp/vasp/Constraining_direction_magnetic_moments.html
1208 : : !!
1209 : : !! SOURCE
1210 : :
1211 : 109 : subroutine mag_penalty(c_dft,mpi_enreg,rhor,nv_constr_dft_r,xred,qgbt,use_gbt)
1212 : :
1213 : : !Arguments ------------------------------------
1214 : : !scalars
1215 : : class(constrained_dft_t),intent(in) :: c_dft
1216 : : integer,intent(in) :: use_gbt
1217 : : real(dp),intent(out) :: nv_constr_dft_r(c_dft%nfftf,c_dft%nspden)
1218 : : type(MPI_type),intent(in) :: mpi_enreg
1219 : : !arrays
1220 : : real(dp),intent(in) :: rhor(c_dft%nfftf,c_dft%nspden)
1221 : : real(dp),intent(in) :: xred(3,c_dft%natom),qgbt(3)
1222 : :
1223 : : !Local variables-------------------------------
1224 : : !scalars
1225 : : integer :: iatom,magconon,natom,nfftf,nspden,ntypat,option
1226 : : integer,parameter :: cplex1=1
1227 : : real(dp):: cmm_x,cmm_y,cmm_z,intgden_proj,norm
1228 : : !arrays
1229 : 109 : real(dp), allocatable :: coeffs_constr_dft(:,:) ! nspden,natom
1230 : 109 : real(dp), allocatable :: intgden(:,:) ! nspden,natom
1231 : 218 : real(dp) :: rhomag(2,c_dft%nspden),spinat_normed(3)
1232 : : ! ***********************************************************************************************
1233 : :
1234 : 109 : magconon=c_dft%magconon
1235 : 109 : natom=c_dft%natom
1236 : 109 : nfftf=c_dft%nfftf
1237 : 109 : nspden=c_dft%nspden
1238 : 109 : ntypat=c_dft%ntypat
1239 : :
1240 [ - + + - : 436 : ABI_MALLOC(coeffs_constr_dft,(nspden,natom))
- + + - +
- + - - +
- + ]
1241 [ - + + - : 327 : ABI_MALLOC(intgden,(nspden,natom))
- + + - -
+ - + ]
1242 : :
1243 : : !We need the integrated magnetic moments and the smoothing function
1244 : : call calcdenmagsph(mpi_enreg,natom,nfftf,c_dft%ngfftf,nspden,ntypat,&
1245 : 109 : c_dft%ratsm,c_dft%ratsph,rhor,c_dft%rprimd,c_dft%typat,xred,1,cplex1,qgbt,use_gbt,intgden=intgden,rhomag=rhomag)
1246 : :
1247 [ + + ]: 218 : call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,[std_out],1,qgbt,c_dft%ratsm,c_dft%ratsph,rhomag,c_dft%typat,c_dft%znucl,c_dft%spinaxis)
1248 : :
1249 : : !Loop over atoms
1250 : : !-------------------------------------------
1251 [ + + ]: 218 : do iatom=1,natom
1252 : :
1253 [ + + ]: 436 : norm = sqrt(sum(c_dft%spinat(:,iatom)**2))
1254 : 109 : spinat_normed(:) = zero
1255 [ + - ]: 109 : if (norm > tol10) then
1256 [ + + ]: 436 : spinat_normed(:) = c_dft%spinat(:,iatom) / norm
1257 [ # # ]: 0 : else if (magconon == 1) then
1258 : : ! if spinat = 0 and we are imposing the direction only, skip this atom
1259 : : cycle
1260 : : end if
1261 : :
1262 : : ! Calculate the x- and y-components of the square bracket term
1263 : 109 : cmm_x = zero
1264 : 109 : cmm_y = zero
1265 : 109 : cmm_z = zero
1266 : 109 : intgden_proj = zero
1267 [ + + ]: 109 : if (nspden == 4) then
1268 [ + + ]: 65 : if (magconon==1) then
1269 : : ! Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector
1270 : : ! This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector
1271 : : intgden_proj=spinat_normed(1)*intgden(2,iatom)+ &
1272 : : & spinat_normed(2)*intgden(3,iatom)+ &
1273 : 54 : & spinat_normed(3)*intgden(4,iatom)
1274 : :
1275 : 54 : cmm_x=intgden(2,iatom)
1276 : 54 : cmm_x=cmm_x-spinat_normed(1)*intgden_proj
1277 : :
1278 : 54 : cmm_y=intgden(3,iatom)
1279 : 54 : cmm_y=cmm_y-spinat_normed(2)*intgden_proj
1280 : :
1281 [ + - ]: 11 : else if (magconon==2 .and. nspden == 4) then
1282 : 11 : cmm_x=intgden(2,iatom)-c_dft%spinat(1,iatom)
1283 : 11 : cmm_y=intgden(3,iatom)-c_dft%spinat(2,iatom)
1284 : : end if
1285 : :
1286 : : ! Calculate the constraining potential for x- and y- components of the mag. mom. vector
1287 : : ! Eric Bousquet has derived the relationship between spin components and potential spin matrix elements:
1288 : : ! 1 = up up = +z
1289 : : ! 2 = down down = -z
1290 : : ! 3 = up down = +x
1291 : : ! 4 = down up = -y
1292 : 65 : coeffs_constr_dft(3,iatom)= 2*c_dft%magcon_lambda*cmm_x
1293 : 65 : coeffs_constr_dft(4,iatom)=-2*c_dft%magcon_lambda*cmm_y
1294 : : end if ! nspden 4
1295 : :
1296 : : ! Calculate the z-component of the square bracket term
1297 [ + + ]: 109 : if (magconon==1) then
1298 [ + + ]: 79 : if (nspden == 4) then
1299 : : ! This apparently enforces the axis of magnetization, not its direction.
1300 : : ! m_z - spinat_z * <m | spinat>
1301 : 54 : cmm_z = intgden(4,iatom) - spinat_normed(3)*intgden_proj
1302 [ + - ]: 25 : else if (nspden == 2) then
1303 : : ! This apparently enforces the direction of the magnetization. So, the behaviour differs in nspden=4 and 2 cases .
1304 : : ! this will be just a sign +/- : are we in the same direction as spinat_z?
1305 : : ! need something more continuous??? To make sure the gradient pushes the state towards FM/AFM?
1306 : 25 : cmm_z = -sign(one, (intgden(1,iatom)-intgden(2,iatom))*spinat_normed(3))
1307 : : end if
1308 [ + - ]: 30 : else if (magconon==2) then
1309 [ + + ]: 30 : if (nspden == 4) then
1310 : 11 : cmm_z=intgden(4,iatom)-c_dft%spinat(3,iatom)
1311 [ + - ]: 19 : else if (nspden == 2) then
1312 : : ! this is up spins - down spins - requested moment ~ 0
1313 : : ! EB: note that intgden comes from calcdenmagsph, which, in nspden=2 case, returns
1314 : : ! intgden(1)=rho_up=n+m
1315 : : ! intgden(2)=rho_dn=n-m
1316 : : ! Then, is the following line be
1317 : : ! cmm_z=half*(intgden(1,iatom)-intgden(2,iatom)) - spinat(3,iatom)
1318 : : ! ??
1319 : 19 : cmm_z=intgden(1,iatom)-intgden(2,iatom) - c_dft%spinat(3,iatom)
1320 : : end if
1321 : : endif
1322 : :
1323 : : ! Calculate the constraining potential for z-component of the mag. mom. vector
1324 : 109 : coeffs_constr_dft(1,iatom)= 2*c_dft%magcon_lambda*cmm_z
1325 : 218 : coeffs_constr_dft(2,iatom)=-2*c_dft%magcon_lambda*cmm_z
1326 : :
1327 : : enddo ! iatom
1328 : :
1329 : : !Now compute the potential in real space
1330 : 109 : option=0
1331 : : call add_atomic_fcts(natom,nspden,c_dft%rprimd,mpi_enreg,nfftf,c_dft%ngfftf,ntypat,option,c_dft%ratsph, &
1332 : 109 : c_dft%ratsm,c_dft%typat,coeffs_constr_dft,nv_constr_dft_r,xred) ! Ben change: add ratsm
1333 : :
1334 [ - + ]: 109 : ABI_FREE(coeffs_constr_dft)
1335 [ - + ]: 109 : ABI_FREE(intgden)
1336 : :
1337 : 109 : end subroutine mag_penalty
1338 : : !!***
1339 : :
1340 : : !!****f* m_dens/mag_penalty_e
1341 : : !! NAME
1342 : : !! mag_penalty_e
1343 : : !!
1344 : : !! FUNCTION
1345 : : !! Compute the energy corresponding to constrained magnetic moments.
1346 : : !!
1347 : : !! INPUTS
1348 : : !! magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size
1349 : : !! spinat=fixed magnetic moments vectors
1350 : : !! magcon_lambda=the size of the penalty terms
1351 : : !!
1352 : : !! OUTPUT
1353 : : !! Epen=penalty contribution to the total energy corresponding to the constrained potential
1354 : : !! Econstr=???
1355 : : !! Eexp=???
1356 : : !!
1357 : : !! SOURCE
1358 : :
1359 : 5 : subroutine mag_penalty_e(magconon,magcon_lambda,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsm,ratsph,rhor,rprimd,spinat,typat,xred,znucl,qgbt,use_gbt,spinaxis)
1360 : :
1361 : : !Arguments ------------------------------------
1362 : : !scalars
1363 : : integer,intent(in) :: natom,magconon,nspden,nfft,ntypat,use_gbt
1364 : : real(dp),intent(in) :: magcon_lambda,ratsm,qgbt(3)
1365 : : !arrays
1366 : : integer, intent(in) :: ngfft(18),typat(natom)
1367 : : real(dp),intent(in) :: spinat(3,natom), rprimd(3,3),znucl(ntypat),spinaxis(3)
1368 : : real(dp),intent(in) :: ratsph(ntypat),rhor(nfft,nspden),xred(3,natom)
1369 : : type(MPI_type),intent(in) :: mpi_enreg
1370 : :
1371 : : !Local variables-------------------------------
1372 : : !scalars
1373 : : integer :: iatom,ii
1374 : : integer,parameter :: cplex1=1 ! dummy argument for calcdenmagsph
1375 : : real(dp) :: intgden_proj, Epen,Econstr,lVp, norm
1376 : : !arrays
1377 : : real(dp) :: intmm(3), mag_1atom(3)
1378 : 5 : real(dp), allocatable :: intgden(:,:)
1379 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),ucvol
1380 : 10 : real(dp) :: rhomag(2,nspden),spinat_normed(3)
1381 : : character(len=500) :: msg
1382 : : ! *********************************************************************
1383 : :
1384 : : !We need the metric because it is needed in calcdenmagsph.F90
1385 : 5 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1386 : :
1387 [ - + + - : 20 : ABI_MALLOC(intgden, (nspden,natom))
- + + - +
- + - - +
- + ]
1388 : :
1389 : : !We need the integrated magnetic moments
1390 : : call calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsm,ratsph,rhor,rprimd,typat,xred,&
1391 : 5 : 1,cplex1,qgbt,use_gbt,intgden=intgden,rhomag=rhomag)
1392 : :
1393 [ + + ]: 10 : call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,[std_out],1,qgbt,ratsm,ratsph,rhomag,typat,znucl,spinaxis)
1394 : :
1395 : 5 : Epen=0
1396 : 5 : Econstr=0
1397 : 5 : lVp=0
1398 : :
1399 : : !Loop over atoms
1400 : : !-------------------------------------------
1401 [ + + ]: 10 : do iatom=1,natom
1402 : :
1403 [ + + ]: 20 : norm = sqrt(sum(spinat(:,iatom)**2))
1404 : 5 : spinat_normed(:) = zero
1405 [ + - ]: 5 : if (norm > tol10) then
1406 [ + + ]: 20 : spinat_normed(:) = spinat(:,iatom) / norm
1407 [ # # ]: 0 : else if (magconon == 1) then
1408 : : ! if spinat = 0 and we are imposing the direction only, skip this atom
1409 : : cycle
1410 : : end if
1411 : : ! Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector
1412 : : ! This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector
1413 : :
1414 : : ! for the collinear spin case, set up a fictitious 3D vector along z
1415 [ + + ]: 5 : if (nspden == 4) then
1416 [ + + ]: 12 : mag_1atom(1:3) = intgden(2:4,iatom)
1417 [ + - ]: 2 : else if (nspden == 2) then
1418 : 2 : mag_1atom = zero
1419 : 2 : mag_1atom(3) = intgden(1,iatom)-intgden(2,iatom)
1420 : : end if
1421 : :
1422 : 5 : intgden_proj = zero
1423 : 5 : intmm = zero
1424 : : ! Calculate the square bracket term
1425 [ + + ]: 5 : if (magconon==1) then
1426 : : intgden_proj=spinat_normed(1)*mag_1atom(1)+ &
1427 : : & spinat_normed(2)*mag_1atom(2)+ &
1428 : 3 : & spinat_normed(3)*mag_1atom(3)
1429 : :
1430 [ + + ]: 12 : do ii=1,3
1431 : 12 : intmm(ii)=mag_1atom(ii)-spinat_normed(ii)*intgden_proj
1432 : : end do
1433 : :
1434 : : ! Calculate the energy Epen corresponding to the constraining potential
1435 : : ! Econstr and lVp do not have a clear meaning (yet)
1436 : 3 : Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3))
1437 : 3 : Econstr=Econstr-magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3))
1438 : 3 : lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3))
1439 : :
1440 [ + - ]: 2 : else if (magconon==2) then
1441 [ + + ]: 8 : do ii=1,3
1442 : 8 : intmm(ii)=mag_1atom(ii)-spinat(ii,iatom)
1443 : : end do
1444 : :
1445 : : ! Calculate the energy Epen corresponding to the constraining potential
1446 : : ! Epen = -Econstr - lVp
1447 : : ! Econstr = -M**2 + spinat**2
1448 : : ! lVp = +2 M \cdot spinat
1449 : 2 : Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3))
1450 : : Econstr=Econstr-magcon_lambda*(mag_1atom(1)*mag_1atom(1)+&
1451 : : & mag_1atom(2)*mag_1atom(2)+&
1452 : : & mag_1atom(3)*mag_1atom(3)) &
1453 : : & +magcon_lambda*(spinat(1,iatom)*spinat(1,iatom)+&
1454 : : & spinat(2,iatom)*spinat(2,iatom)+&
1455 : 2 : & spinat(3,iatom)*spinat(3,iatom))
1456 : 2 : lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3))
1457 : : end if
1458 : :
1459 : 5 : write(msg, *) 'atom constraining magnetic field'
1460 : 5 : call wrtout(std_out,msg)
1461 : : write(msg, '(I3,A2,E12.5,A2,E12.5,A2,E12.5)') &
1462 : 5 : iatom,' ',magcon_lambda*intmm(1),' ',magcon_lambda*intmm(2),' ',magcon_lambda*intmm(3)
1463 : 10 : call wrtout(std_out,msg)
1464 : :
1465 : : ! End loop over atoms
1466 : : ! -------------------------------------------
1467 : : end do
1468 : :
1469 : : !Printing
1470 : 5 : write(msg, '(A17,E10.3)' ) ' magcon_lambda = ',magcon_lambda
1471 : 5 : call wrtout(std_out,msg)
1472 : 5 : write(msg, '(A17,E12.5)' ) ' Lagrange penalty = ',Epen
1473 : 5 : call wrtout(std_out,msg)
1474 : 5 : write(msg, '(A17,E12.5)' ) ' E_constraint = ',Econstr
1475 : 5 : call wrtout(std_out,msg)
1476 : 5 : write(msg, '(A17,E12.5)' ) ' lVp = ',lVp
1477 : 5 : call wrtout(std_out,msg)
1478 : :
1479 [ - + ]: 5 : ABI_FREE(intgden)
1480 : :
1481 : 5 : end subroutine mag_penalty_e
1482 : : !!***
1483 : :
1484 : : !!****f* m_dens/calcdenmagsph
1485 : : !! NAME
1486 : : !! calcdenmagsph
1487 : : !!
1488 : : !! FUNCTION
1489 : : !! Compute and print integral of total density inside spheres around atoms,
1490 : : !! or optionally of integral of potential residual.
1491 : : !! Also can compute the contributions to forces and stresses due to density-magnetization type constraints.
1492 : : !!
1493 : : !! INPUTS
1494 : : !! mpi_enreg=information about MPI parallelization
1495 : : !! natom=number of atoms in cell.
1496 : : !! nfft=(effective) number of FFT grid points (for this processor)
1497 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1498 : : !! nspden=number of spin-density components
1499 : : !! ntypat=number of atom types
1500 : : !! ratopt= if 1 the atomic spheres are defined in real space
1501 : : !! if 2 the atomic spheres are dfined in reciprocal space and then Fourier transformed
1502 : : !! option = if not larger than 10, then a density is input , if larger than 10 then a potential residual is input.
1503 : : !! ratsm=smearing width for ratsph
1504 : : !! ratsph(ntypat)=radius of spheres around atoms
1505 : : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
1506 : : !! (total in first half and spin-up in second half if nspden=2)
1507 : : !! (total in first comp. and magnetization in comp. 2 to 4 if nspden=4)
1508 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1509 : : !! typat(natom)=type of each atom
1510 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
1511 : : !! [qphon(3)]= perturbation wave vector.
1512 : : !!
1513 : : !! OUTPUT
1514 : : !! dentot(nspden)=integrated density (magnetization...) over full u.c. vol. Optional argument
1515 : : !! gr_intgden(3,nspden,natom)=grad wrt atomic positions, of integrated density (magnetization...) for each atom in a sphere. Optional arg
1516 : : !! intgden(nspden, natom)=integrated density (magnetization...) for each atom in a sphere of radius ratsph. Optional arg
1517 : : !! Note that when intgden is present, the definition of the spherical integration function changes, as it is smoothed.
1518 : : !! intgf2(natom,natom)=overlaps of the spherical integration functions for each atom in a sphere of radius ratsph. Optional arg
1519 : : !! rhomag(2,nspden)=integrated complex density (magnetization...) over full u.c. vol. Optional argument
1520 : : !! In collinear case component 1 is total density and 2 is _magnetization_ up-down
1521 : : !! In non collinear case component 1 is total density, and 2:4 are the magnetization vector
1522 : : !! strs_intgden(6,nspden,natom)=stress contribution due to constrained integrated density (magnetization...), due to each atom. Optional arg
1523 : : !! fatsph(nfft,natom)= functions defining the atomic spheres of integration in real space
1524 : : !! taumr(nfft,natom,3)= array describing r-xred(iatom) at any point of the FFT grid
1525 : : !! Rest is printing
1526 : : !!
1527 : : !! SOURCE
1528 : :
1529 : 57084 : subroutine calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsm,ratsph,rhor,rprimd,typat,xred,&
1530 : 21740 : & option,cplex,qgbt,use_gbt,dentot,gr_intgden,intgden,intgf2,rhomag,strs_intgden,fatsph,qphon,taumr)
1531 : :
1532 : : !Arguments ---------------------------------------------
1533 : : !scalars
1534 : : integer,intent(in) :: natom,nfft,nspden,ntypat
1535 : : real(dp),intent(in) :: ratsm
1536 : : type(MPI_type),intent(in) :: mpi_enreg
1537 : : integer,intent(in) :: option
1538 : : integer,intent(in) :: cplex
1539 : : integer,intent(in) :: use_gbt
1540 : : !arrays
1541 : : integer,intent(in) :: ngfft(18),typat(natom)
1542 : : real(dp),intent(in) :: ratsph(ntypat),rhor(cplex*nfft,nspden),rprimd(3,3)
1543 : : real(dp),intent(in) :: xred(3,natom)
1544 : : real(dp),intent(in) :: qgbt(3)
1545 : : real(dp),intent(out),optional :: dentot(nspden)
1546 : : real(dp),intent(out),optional :: gr_intgden(3,nspden,natom)
1547 : : real(dp),intent(out),optional :: intgden(cplex,nspden,natom)
1548 : : real(dp),intent(out),optional :: intgf2(natom,natom)
1549 : : real(dp),intent(out),optional :: rhomag(2,nspden)
1550 : : real(dp),intent(out),optional :: strs_intgden(6,nspden,natom)
1551 : : real(dp),intent(out),optional,target :: fatsph(nfft,natom)
1552 : : real(dp),intent(in),optional :: qphon(3)
1553 : : real(dp),intent(out),optional,target :: taumr(nfft,natom,3)
1554 : : !Local variables ------------------------------
1555 : :
1556 : : !scalars
1557 : : integer,parameter :: ndir=3,ishift=5
1558 : : integer :: i1,i2,i3,iatom,ierr,ifft_local,ii,isp,ispden,ix,iy,iz,izloc,jatom,n1,n1a,n1b,n2,ifft,ifft_local_cplex
1559 : : integer :: neighbor_overlap,n2a,n2b,n3,n3a,n3b,nfftot,n4,n5,n6
1560 : : ! integer :: n1c, n2c, n3c
1561 : : integer :: jfft
1562 : : real(dp) :: arg,phr1d_im,phr1d_re
1563 : : real(dp),parameter :: delta=0.99_dp
1564 : 38056 : real(dp) :: difx,dify,difz,r2,r2atsph,rr1,rr2,rr3,rx,ry,rz,qr,mx,my,mz,rhor_local(nspden)
1565 : : real(dp) :: dfsm,fact,fsm,ratsm2,ucvol
1566 : : logical :: grid_found
1567 : : !arrays
1568 : 19028 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1569 : 38056 : integer :: overlap_ij(natom,natom)
1570 : : real(dp) :: gmet(3,3),gprimd(3,3),gr_intg(3,4)
1571 : 38056 : real(dp) :: intg(cplex,4),qphon_(3),rhomag_(2,nspden)
1572 : : ! real(dp) :: intg_im(4),intg_re(4)
1573 : : real(dp) :: strs(3,3),strs_cartred(3,3),strs_intg(6,4),tsec(2)
1574 : 57084 : real(dp) :: dist_ij(natom,natom),intgden_(cplex,nspden,natom)!,intgden_im_(nspden,natom)
1575 : 38056 : real(dp) :: my_xred(3, natom), rmet(3,3),xshift(3, natom), taumr_local(3)
1576 : 19028 : real(dp), allocatable :: fsm_atom(:,:)
1577 : 19028 : real(dp), contiguous, pointer :: fatsph_(:,:),taumr_(:,:,:)
1578 : : !real(dp) :: rprimd_mod(3,3),strain
1579 : : ! *************************************************************************
1580 : :
1581 : : !MG NOTE: the computation of intg is clearly wrong when cplex = 2 (DFPT)
1582 : :
1583 : 19028 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
1584 : 19028 : n4=ngfft(4);n5=ngfft(5);n6=ngfft(6)
1585 : 19028 : nfftot=n1*n2*n3
1586 : :
1587 : : !Manage optinal arguments
1588 [ + + + + : 261027 : if(present(intgden)) intgden=zero
+ + ]
1589 [ + + + + : 65502 : if(present(gr_intgden)) gr_intgden=zero
+ + + + ]
1590 [ + + + + : 98016 : if(present(strs_intgden)) strs_intgden=zero
+ + + + ]
1591 : 19028 : qphon_=zero
1592 [ + + ]: 19028 : if(present(qphon))then
1593 : 74 : qphon_=qphon
1594 : : endif
1595 : :
1596 [ + + ]: 19028 : if(present(fatsph)) then
1597 : 62 : fatsph_ => fatsph
1598 [ + + + + ]: 952754 : fatsph_=zero
1599 : : end if
1600 [ + + ]: 19028 : if(present(taumr)) then
1601 : 62 : taumr_ => taumr
1602 [ + + + + : 2858324 : taumr_=zero
+ + ]
1603 : : end if
1604 : :
1605 : 19028 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1606 : :
1607 : : ! This routine is not able to handle xred positions that are "far" from the
1608 : : ! first unit cell so wrap xred into [0, 1[ interval here.
1609 [ + + + + ]: 204944 : call wrap2_zero_one(xred, my_xred, xshift)
1610 : :
1611 : : !If intgf2 present, check that the spheres do not overlap. If they overlap, one needs to treat explicitly the neighbors.
1612 [ + + ]: 19028 : if(present(intgf2))then
1613 : 180 : neighbor_overlap=0
1614 [ + + + + ]: 498 : dist_ij(:,:)=zero
1615 [ + + + + ]: 498 : intgf2(:,:)=zero
1616 [ + + + + ]: 498 : overlap_ij(:,:)=0
1617 : 54 : dist_ij(1,1)=dist2(xred(:,1),xred(:,1),rprimd,-1)
1618 [ + + ]: 180 : do iatom=1,natom
1619 : 126 : dist_ij(iatom,iatom)=dist_ij(1,1)
1620 : 126 : overlap_ij(iatom,iatom)=1
1621 [ + + ]: 402 : do jatom=iatom,natom
1622 [ + + ]: 222 : if(iatom/=jatom)dist_ij(iatom,jatom)=dist2(xred(:,iatom),xred(:,jatom),rprimd,1)
1623 [ + + ]: 348 : if(dist_ij(iatom,jatom)-ratsph(typat(iatom))-ratsph(typat(jatom))<tol10)then
1624 : 49 : overlap_ij(iatom,jatom)=1
1625 : 49 : neighbor_overlap=1
1626 : : endif
1627 : : enddo
1628 : : enddo
1629 [ + + ]: 54 : if(neighbor_overlap==1)then
1630 [ - + + - : 76 : ABI_MALLOC(fsm_atom,(nfft,natom))
- + + - +
- + - - +
- + ]
1631 [ + + + + ]: 626181 : fsm_atom(:,:)=zero
1632 : : endif
1633 : : endif
1634 : :
1635 : : !Get the distrib associated with this fft_grid
1636 : 19028 : grid_found=.false.
1637 : :
1638 [ + + ]: 19028 : if(n2 == mpi_enreg%distribfft%n2_coarse ) then
1639 [ + + ]: 13669 : if(n3== size(mpi_enreg%distribfft%tab_fftdp3_distrib)) then
1640 : 13615 : fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3_distrib
1641 : 13615 : ffti3_local => mpi_enreg%distribfft%tab_fftdp3_local
1642 : 13615 : grid_found=.true.
1643 : : end if
1644 : : end if
1645 : :
1646 [ + + ]: 19028 : if(n2 == mpi_enreg%distribfft%n2_fine ) then
1647 [ + - ]: 5621 : if(n3 == size(mpi_enreg%distribfft%tab_fftdp3dg_distrib)) then
1648 : 5621 : fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3dg_distrib
1649 : 5621 : ffti3_local => mpi_enreg%distribfft%tab_fftdp3dg_local
1650 : : grid_found = .true.
1651 : : end if
1652 : : end if
1653 : :
1654 [ + - ]: 13407 : ABI_CHECK(grid_found, "Unable to find an allocated distrib for this fft grid")
1655 : :
1656 : : !Loop over atoms
1657 : : !-------------------------------------------
1658 [ + + ]: 65507 : do iatom=1,natom
1659 : :
1660 : : ! Define a "box" around the atom that extends outside of the unit cell
1661 : 46479 : r2atsph=1.0000001_dp*ratsph(typat(iatom))**2
1662 : 46479 : rr1=sqrt(r2atsph*gmet(1,1))
1663 : 46479 : rr2=sqrt(r2atsph*gmet(2,2))
1664 : 46479 : rr3=sqrt(r2atsph*gmet(3,3))
1665 : :
1666 : 46479 : n1a=int((my_xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1
1667 : 46479 : n1b=int((my_xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1
1668 : 46479 : n2a=int((my_xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2
1669 : 46479 : n2b=int((my_xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2
1670 : 46479 : n3a=int((my_xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3
1671 : 46479 : n3b=int((my_xred(3,iatom)+rr3+ishift)*n3 )-ishift*n3
1672 : :
1673 : : !This is the "width" of the zone of smearing, in term of the square of radius
1674 : 46479 : ratsm2 = (2*ratsph(typat(iatom))-ratsm)*ratsm
1675 : :
1676 [ + + + + ]: 421259 : intg(:,:)=zero
1677 : 46479 : gr_intg(:,:)=zero
1678 : 46479 : strs_intg(:,:)=zero
1679 : :
1680 [ + + ]: 644300 : do i3=n3a,n3b
1681 : 597821 : iz=mod(i3+ishift*n3,n3)
1682 : :
1683 [ + + ]: 644300 : if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
1684 : :
1685 : 586685 : izloc = ffti3_local(iz+1) - 1
1686 : 586685 : difz=dble(i3)/dble(n3)-my_xred(3,iatom)
1687 [ + + ]: 10439580 : do i2=n2a,n2b
1688 : 9852895 : iy=mod(i2+ishift*n2,n2)
1689 : 9852895 : dify=dble(i2)/dble(n2)-my_xred(2,iatom)
1690 [ + + ]: 231492731 : do i1=n1a,n1b
1691 : 221053151 : ix=mod(i1+ishift*n1,n1)
1692 : 221053151 : difx=dble(i1)/dble(n1)-my_xred(1,iatom)
1693 : : !DEBUG
1694 : : ! if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
1695 : : ! difx=dble(i1)/dble(n1)-(my_xred(1,iatom)+0.00005)
1696 : : ! endif
1697 : : !ENDDEBUG
1698 : 221053151 : rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
1699 : :
1700 : : !DEBUG
1701 : : ! if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
1702 : : ! strain=-0.001
1703 : : ! rx=difx*rprimd(1,1)*(one+strain)+dify*rprimd(1,2)+difz*rprimd(1,3)
1704 : : ! endif
1705 : : !ENDDEBUG
1706 : :
1707 : 221053151 : ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1708 : 221053151 : rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1709 : 221053151 : r2=rx**2+ry**2+rz**2
1710 : :
1711 : : ! Identify the fft indexes of the rectangular grid around the atom
1712 [ + + ]: 221053151 : if(r2 > r2atsph) then
1713 : : cycle
1714 : : end if
1715 : :
1716 : 93898393 : call radsmear(dfsm,fsm,r2,r2atsph,ratsm2)
1717 : :
1718 : 93898393 : ifft_local=1+ix+n1*(iy+n2*izloc)
1719 : :
1720 [ + + ]: 93898393 : if(present(fatsph)) then
1721 : 60636 : fatsph_(ifft_local,iatom)=fsm
1722 : : end if
1723 [ + + ]: 93898393 : if(present(taumr)) then
1724 : 60636 : taumr_(ifft_local,iatom,1)=difx
1725 : 60636 : taumr_(ifft_local,iatom,2)=dify
1726 : 60636 : taumr_(ifft_local,iatom,3)=difz
1727 : : end if
1728 : :
1729 : : ! Compute the finite-q real-space phase
1730 : 93898393 : taumr_local(1)=difx
1731 : 93898393 : taumr_local(2)=dify
1732 : 93898393 : taumr_local(3)=difz
1733 [ + + ]: 375593572 : arg=two_pi*dot_product(qphon_,taumr_local)
1734 : 93898393 : phr1d_re=dcos(arg)
1735 : 93898393 : phr1d_im=dsin(arg)
1736 : 93898393 : ifft_local_cplex=1+cplex*(ifft_local-1)
1737 : :
1738 [ + + ]: 93898393 : if(present(intgf2))then
1739 : : ! intgden_(1,iatom)= integral of the square of the spherical integrating function
1740 [ + + ]: 109724 : if(neighbor_overlap==0)intgf2(iatom,iatom)=intgf2(iatom,iatom)+fsm*fsm
1741 : 13136 : if(neighbor_overlap==1)fsm_atom(ifft_local,iatom)=fsm_atom(ifft_local,iatom)+fsm
1742 : : endif
1743 : :
1744 [ + + + + ]: 93898393 : if (nspden==4 .and. use_gbt /= 0) then
1745 : : ! If GBT is on, we have to include the e^{iq.r} phase
1746 : 750873 : qr = two_pi * (qgbt(1)*dble(i1)/dble(n1) + qgbt(2)*dble(i2)/dble(n2) + qgbt(3)*dble(i3)/dble(n3))
1747 : 750873 : rhor_local(1) = rhor(ifft_local,1)
1748 : 750873 : mx = rhor(ifft_local,2)
1749 : 750873 : my = rhor(ifft_local,3)
1750 : 750873 : mz = rhor(ifft_local,4)
1751 : : ! e^-iqr :cos(qr)-isin(qr)
1752 : 750873 : rhor_local(2) = cos(qr)*mx - sin(qr)*my ! mx
1753 : 750873 : rhor_local(3) = sin(qr)*mx + cos(qr)*my ! my
1754 : 750873 : rhor_local(4) = mz
1755 : : else
1756 [ + + ]: 272759490 : rhor_local(1:nspden) = rhor(ifft_local,1:nspden)
1757 : : end if
1758 : : ! Integral of density or potential residual
1759 [ + + ]: 93898393 : if (cplex==1) then
1760 : : ! intg(1,1:nspden)=intg(1,1:nspden)+fsm*rhor(ifft_local,1:nspden)
1761 [ + + ]: 274552096 : intg(1,1:nspden)=intg(1,1:nspden)+fsm*rhor_local(1:nspden)
1762 [ + - ]: 457533 : else if (cplex==2) then
1763 [ + + + + ]: 1830132 : if (sum(qphon_(:)**2)<tol8) then
1764 [ + + ]: 1823959 : intg(1,1:nspden)=intg(1,1:nspden)+fsm*rhor(2*ifft_local-1,1:nspden)
1765 [ + + ]: 1823959 : intg(2,1:nspden)=intg(2,1:nspden)+fsm*rhor(2*ifft_local ,1:nspden)
1766 : : else
1767 [ + + ]: 137800 : intg(1,1:nspden)=intg(1,1:nspden)+phr1d_re*fsm*rhor(2*ifft_local-1,1:nspden)-phr1d_im*fsm*rhor(2*ifft_local ,1:nspden)
1768 [ + + ]: 137800 : intg(2,1:nspden)=intg(2,1:nspden)+phr1d_re*fsm*rhor(2*ifft_local ,1:nspden)+phr1d_im*fsm*rhor(2*ifft_local-1,1:nspden)
1769 : : end if
1770 : : end if
1771 [ + + + - : 197649681 : if((present(gr_intgden).or.present(strs_intgden)).and. option<10 .and. ratsm2>tol12)then
+ - ]
1772 [ + + ]: 10676076 : do ispden=1,nspden
1773 : : ! For cplex=1, rhor_local is the GBT phase-rotated density
1774 : : ! (and equals rhor when GBT is off).
1775 : : ! The GBT phase rotation is not implemented for cplex=2
1776 [ + - ]: 7343506 : if (cplex == 1) then
1777 : 7343506 : fact=dfsm*rhor_local(ispden)
1778 : : else
1779 : 0 : fact=dfsm*rhor(ifft_local_cplex,ispden)
1780 : : end if
1781 : : ! fact=dfsm*rhor(ifft_local_cplex,ispden)
1782 [ + - ]: 7343506 : if(present(gr_intgden))then
1783 : 7343506 : gr_intg(1,ispden)=gr_intg(1,ispden)+difx*fact
1784 : 7343506 : gr_intg(2,ispden)=gr_intg(2,ispden)+dify*fact
1785 : 7343506 : gr_intg(3,ispden)=gr_intg(3,ispden)+difz*fact
1786 : : endif
1787 [ + - ]: 10676076 : if(present(strs_intgden))then
1788 : 7343506 : strs_intg(1,ispden)=strs_intg(1,ispden)+difx*difx*fact
1789 : 7343506 : strs_intg(2,ispden)=strs_intg(2,ispden)+dify*dify*fact
1790 : 7343506 : strs_intg(3,ispden)=strs_intg(3,ispden)+difz*difz*fact
1791 : 7343506 : strs_intg(4,ispden)=strs_intg(4,ispden)+dify*difz*fact
1792 : 7343506 : strs_intg(5,ispden)=strs_intg(5,ispden)+difx*difz*fact
1793 : 7343506 : strs_intg(6,ispden)=strs_intg(6,ispden)+difx*dify*fact
1794 : : endif
1795 : : enddo
1796 : : endif
1797 : : end do
1798 : : end do
1799 : : end if
1800 : : end do
1801 : :
1802 : :
1803 : : !DEBUG
1804 : : ! n1c=(n1b+n1a)/2
1805 : : ! n2c=(n2b+n2a)/2
1806 : : ! do i3= n3a-5,n3b+5
1807 : : ! n1c=mod(n1c+ishift*n1,n1)
1808 : : ! n2c=mod(n2c+ishift*n2,n2)
1809 : : ! iz=mod(i3+ishift*n3,n3)
1810 : : ! end do
1811 : : !ENDDEBUG
1812 : :
1813 : :
1814 [ + + ]: 46479 : if(present(intgf2) .and. neighbor_overlap==0)then
1815 : 76 : intgf2(iatom,iatom)=intgf2(iatom,iatom)*ucvol/dble(nfftot)
1816 : : endif
1817 : :
1818 [ + + + + ]: 421259 : intg(:,:)=intg(:,:)*ucvol/dble(nfftot)
1819 : :
1820 [ + + + - : 46479 : if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
+ - ]
1821 : : ! Convert to gradient in reduced coordinates
1822 [ + + + + : 165466 : gr_intg=matmul(rmet,gr_intg)
+ + ]
1823 [ + + + + ]: 53074 : gr_intg(:,:)=-gr_intg(:,:)*two*ucvol/dble(nfftot)
1824 : : !DEBUG
1825 : : ! write(6,*)' calcdenmagsph : intg(1)=',intg(1)
1826 : : ! write(6,*)' calcdenmagsph : iatom,gr_intg(1,1)=',iatom,gr_intg(1,1)
1827 : : ! call flush(6)
1828 : : ! stop
1829 : : !ENDDEBUG
1830 : : endif
1831 : :
1832 [ + + + - : 46479 : if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then
+ - ]
1833 : : ! Convert to stress in cartesian coordinates, for each spin constraint
1834 [ + + ]: 15610 : do isp=1,4
1835 : : ! First change the representation
1836 : 12488 : strs(1,1)=strs_intg(1,isp) ; strs(2,2)=strs_intg(2,isp) ; strs(3,3)=strs_intg(3,isp)
1837 : 12488 : strs(2,3)=strs_intg(4,isp) ; strs(3,2)=strs_intg(4,isp)
1838 : 12488 : strs(1,3)=strs_intg(5,isp) ; strs(3,1)=strs_intg(5,isp)
1839 : 12488 : strs(1,2)=strs_intg(6,isp) ; strs(2,1)=strs_intg(6,isp)
1840 : :
1841 : : ! Then perform representation change, following Eq.(25) in Hamann2005
1842 [ + + ]: 49952 : do ii=1,3
1843 [ + + + + : 611912 : strs_cartred(:,ii)=matmul(rprimd,strs(:,ii))
+ + ]
1844 : : enddo
1845 [ + + ]: 49952 : do ii=1,3
1846 [ + + + + : 611912 : strs(ii,:)=matmul(rprimd,strs_cartred(ii,:))
+ + ]
1847 : : enddo
1848 : :
1849 : : !DEBUG
1850 : : ! rprimd_mod=rprimd
1851 : : ! rprimd_mod(1,1)=rprimd(1,1)*(one+strain)
1852 : : ! do ii=1,3
1853 : : ! strs_cartred(:,ii)=matmul(rprimd_mod,strs(:,ii))
1854 : : ! enddo
1855 : : ! do ii=1,3
1856 : : ! strs(ii,:)=matmul(rprimd_mod,strs_cartred(ii,:))
1857 : : ! enddo
1858 : : !ENDDEBUG
1859 : :
1860 : 12488 : strs_intg(1,isp)=two*strs(1,1) ; strs_intg(2,isp)=two*strs(2,2) ; strs_intg(3,isp)=two*strs(3,3)
1861 : 12488 : strs_intg(4,isp)=strs(2,3)+strs(3,2)
1862 : 12488 : strs_intg(5,isp)=strs(1,3)+strs(3,1)
1863 : 15610 : strs_intg(6,isp)=strs(1,2)+strs(2,1)
1864 : : enddo
1865 [ + + + + ]: 90538 : strs_intg(:,:)=strs_intg(:,:)/dble(nfftot)
1866 : : !DEBUG
1867 : : ! strs_intg(:,:)=zero
1868 : : ! write(6,*)' calcdenmagsph : iatom,-strs_intg(1,1)*ucvol=',iatom,strs_intg(1,1)*ucvol
1869 : : ! call flush(6)
1870 : : ! stop
1871 : : !ENDDEBUG
1872 : : endif
1873 : :
1874 [ + + + + ]: 65507 : if(nspden==2 .and. option/=11)then
1875 : : ! Specific treatment of collinear density, due to the storage mode.
1876 : : ! intgden_(1,iatom)= integral of up density
1877 : : ! intgden_(2,iatom)= integral of dn density
1878 [ + + ]: 30419 : intgden_(1:cplex,1,iatom)=intg(1:cplex,2)
1879 [ + + ]: 30419 : intgden_(1:cplex,2,iatom)=intg(1:cplex,1)-intg(1:cplex,2)
1880 : : !if (cplex==2) then
1881 : : ! intgden_im_(1,iatom)=intg(2,2)
1882 : : ! intgden_im_(2,iatom)=intg(2,1)-intg(2,2)
1883 : : !endif
1884 [ + + + - : 15149 : if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
+ - ]
1885 [ + + ]: 2040 : gr_intgden(:,1,iatom)=gr_intg(:,2)
1886 [ + + ]: 2040 : gr_intgden(:,2,iatom)=gr_intg(:,1)-gr_intg(:,2)
1887 : : endif
1888 [ + + + - : 15149 : if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then
+ - ]
1889 [ + + ]: 3570 : strs_intgden(:,1,iatom)=strs_intg(:,2)
1890 [ + + ]: 3570 : strs_intgden(:,2,iatom)=strs_intg(:,1)-strs_intg(:,2)
1891 : : endif
1892 : : else
1893 [ + + + + ]: 167044 : intgden_(1:cplex,1:nspden,iatom)=intg(1:cplex,1:nspden)
1894 [ + + + - : 31330 : if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
+ - ]
1895 [ + + + + ]: 41884 : gr_intgden(:,1:nspden,iatom)=gr_intg(:,1:nspden)
1896 : : endif
1897 [ + + + - : 31330 : if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then
+ - ]
1898 [ + + + + ]: 71338 : strs_intgden(:,1:nspden,iatom)=strs_intg(:,1:nspden)
1899 : : endif
1900 : : endif
1901 : :
1902 : : end do ! iatom
1903 : :
1904 : : !-------------------------------------------
1905 : : !
1906 : : ! In case intgf2 must be computed, while the atoms overlap, a double loop over atoms is needed
1907 [ + + ]: 19028 : if(present(intgf2) .and. neighbor_overlap==1)then
1908 [ + + ]: 69 : do iatom=1,natom
1909 [ + + ]: 149 : do jatom=iatom,natom
1910 [ + - ]: 149 : if(overlap_ij(iatom,jatom)/=0)then
1911 [ + + ]: 1935 : do i3=1,n3
1912 : 1836 : iz=mod(i3,n3)
1913 [ + - ]: 1935 : if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
1914 : 1836 : izloc = ffti3_local(iz+1) - 1
1915 [ + + ]: 42300 : do i2=1,n2
1916 : 40464 : iy=mod(i2,n2)
1917 [ + + ]: 1128348 : do i1=1,n1
1918 : 1086048 : ix=mod(i1,n1)
1919 : 1086048 : ifft_local=1+ix+n1*(iy+n2*izloc)
1920 : 1126512 : intgf2(iatom,jatom)=intgf2(iatom,jatom)+fsm_atom(ifft_local,iatom)*fsm_atom(ifft_local,jatom)
1921 : : enddo
1922 : : enddo
1923 : : endif
1924 : : enddo ! i3
1925 : 99 : intgf2(iatom,jatom)=intgf2(iatom,jatom)*ucvol/dble(nfftot)
1926 : : endif
1927 : : enddo
1928 [ + + ]: 69 : if(iatom/=1)then
1929 [ + + ]: 80 : do jatom=1,iatom-1
1930 : 80 : intgf2(iatom,jatom)=intgf2(jatom,iatom)
1931 : : enddo
1932 : : endif
1933 : : enddo
1934 [ - + ]: 19 : ABI_FREE(fsm_atom)
1935 : : endif
1936 : :
1937 : : !MPI parallelization
1938 [ + + ]: 19028 : if(present(intgf2)) then
1939 [ - + ]: 54 : if(mpi_enreg%nproc_fft>1)then
1940 : 0 : call timab(48,1,tsec)
1941 : 0 : call xmpi_sum(intgf2,mpi_enreg%comm_fft,ierr)
1942 : 0 : call timab(48,2,tsec)
1943 : : end if
1944 : : end if
1945 : :
1946 : : !-------------------------------------------
1947 : :
1948 : : !MPI parallelization
1949 [ + + - + ]: 19028 : if(present(intgden) .or. option/=0) then
1950 [ + + ]: 18974 : if(mpi_enreg%nproc_fft>1)then
1951 : 556 : call timab(48,1,tsec)
1952 : 556 : call xmpi_sum(intgden_,mpi_enreg%comm_fft,ierr)
1953 : 556 : call timab(48,2,tsec)
1954 : : end if
1955 [ + - + + : 261027 : if(present(intgden))intgden = intgden_
+ + + + ]
1956 : : end if
1957 : :
1958 [ + + + - : 19028 : if(present(gr_intgden) .and. option<10 .and. ratsm2>tol12) then
+ - ]
1959 [ - + ]: 1294 : if(mpi_enreg%nproc_fft>1)then
1960 : 0 : call timab(48,1,tsec)
1961 : 0 : call xmpi_sum(gr_intgden,mpi_enreg%comm_fft,ierr)
1962 : 0 : call timab(48,2,tsec)
1963 : : end if
1964 : : end if
1965 : :
1966 [ + + + - : 19028 : if(present(strs_intgden) .and. option<10 .and. ratsm2>tol12) then
+ - ]
1967 [ - + ]: 1294 : if(mpi_enreg%nproc_fft>1)then
1968 : 0 : call timab(48,1,tsec)
1969 : 0 : call xmpi_sum(strs_intgden,mpi_enreg%comm_fft,ierr)
1970 : 0 : call timab(48,2,tsec)
1971 : : end if
1972 : : end if
1973 : :
1974 : : !EB - Compute magnetization of the whole cell
1975 [ + + ]: 19028 : if(present(dentot) .or. present(rhomag))then
1976 [ + + + + ]: 74359 : rhomag_(:,:)=zero
1977 [ + + ]: 11341 : if(nspden==2) then
1978 [ + + ]: 54959314 : do ifft=1,nfft
1979 : 54956654 : jfft=1+cplex*(ifft-1)
1980 [ + + ]: 110314148 : rhomag_(1:cplex,1)=rhomag_(1:cplex,1)+rhor(jfft:jfft+cplex-1,1) ! real & imag part of density
1981 [ + + ]: 110316808 : rhomag_(1:cplex,2)=rhomag_(1:cplex,2)+2*rhor(jfft:jfft+cplex-1,2)-rhor(jfft:jfft+cplex-1,1) ! real & imag part of magnetization
1982 : : end do
1983 [ + + ]: 8681 : else if(nspden==4) then
1984 [ + + ]: 27301758 : do ifft=1,nfft
1985 : 27299423 : jfft=1+cplex*(ifft-1)
1986 [ + + + + ]: 248593014 : rhomag_(1:cplex,1:nspden)=rhomag_(1:cplex,1:nspden)+rhor(jfft:jfft+cplex-1,1:nspden)
1987 : : end do
1988 : : end if
1989 : :
1990 [ + + + + ]: 54071 : rhomag_(1:cplex,1:nspden)=rhomag_(1:cplex,1:nspden)*ucvol/dble(nfftot)
1991 : :
1992 : : !MPI parallelization
1993 [ + + ]: 11341 : if(mpi_enreg%nproc_fft>1)then
1994 : 444 : call timab(48,1,tsec)
1995 : 444 : call xmpi_sum(rhomag_,mpi_enreg%comm_fft,ierr)
1996 : 444 : call timab(48,2,tsec)
1997 : : end if
1998 : :
1999 [ + + + + ]: 14951 : if(present(dentot)) dentot(:)=rhomag_(1,:)
2000 [ + - + + : 74359 : if(present(rhomag)) rhomag(:,:)=rhomag_(:,:)
+ + ]
2001 : : endif
2002 : :
2003 : : !DEBUG BUT KEEP
2004 : : !if(.false.) call printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,'DEN.vtk')
2005 : :
2006 [ + + + + : 59796 : end subroutine calcdenmagsph
+ + + + +
+ + + + +
- + ]
2007 : : !!***
2008 : :
2009 : : !!****f* m_dens/prtdenmagsph
2010 : : !! NAME
2011 : : !! prtdenmagsph
2012 : : !!
2013 : : !! FUNCTION
2014 : : !! Print integral of total density inside spheres around atoms,
2015 : : !! and optionally integral of potential residual (also gradient of the constraint energy wrt constraint).
2016 : : !!
2017 : : !! INPUTS
2018 : : !! intgden(nspden, natom)=integrated rhor or potential residual, for each atom in a sphere of radius ratsph.
2019 : : !! if option <10, intgden is a density (+magnetization)
2020 : : !! Representation differs according to nspden :
2021 : : !! if nspden=1, total density
2022 : : !! if nspden=2, spin up, then spin down
2023 : : !! if nspden=4, total density, then mag_x, mag_y, mag_z
2024 : : !! if 20>option>=10, intgden is a potential residual (+spin magnetic field residual).
2025 : : !! Representation is: first, mean potential; then (if nspden>=2) B_z for nspden=2, B_x, B_y and B_z for nspden=4.
2026 : : !! if option>=20, intgden is a a gradient wrt target (=torque), also potential residual (+spin magnetic field residual) multiplied by f-1.
2027 : : !! natom=number of atoms in cell.
2028 : : !! nspden=number of spin-density components
2029 : : !! ntypat=number of atom types
2030 : : !! units=unit numbers for printing
2031 : : !! option = if not larger than 10, then a density is input , if between 10 and 19 then a potential residual is input, if beyond, a torque is input..
2032 : : !! When 1, 11, 21, the default printing is on (to unit nunit), if -1, 2, 3, 4, special printing options, if 0 no printing.
2033 : : !! ratsm=smearing width for ratsph
2034 : : !! ratsph(ntypat)=radius of spheres around atoms
2035 : : !! rhomag(2,nspden)=integral of charge or magnetization over the whole cell (also taking into account a possible imaginary part for DFPT).
2036 : : !! typat(natom)=type of each atom
2037 : : !! ziontypat(ntypat)= --optional-- ionic charge of each atomic type
2038 : : !!
2039 : : !! OUTPUT
2040 : : !! Printing
2041 : : !!
2042 : : !! SOURCE
2043 : :
2044 [ - + ]: 9823 : subroutine prtdenmagsph(cplex, intgden, natom, nspden, ntypat, units, option, qgbt, ratsm, ratsph, rhomag, typat, znucl, spinaxis, ziontypat)
2045 : :
2046 : : !Arguments ---------------------------------------------
2047 : : !scalars
2048 : : integer,intent(in) :: natom,nspden,ntypat
2049 : : integer,intent(in) :: units(:)
2050 : : real(dp),intent(in) :: ratsm
2051 : : integer ,intent(in) :: option
2052 : : integer, intent(in) :: cplex
2053 : : !arrays
2054 : : integer,intent(in) :: typat(natom)
2055 : : real(dp),intent(in) :: intgden(cplex,nspden,natom),qgbt(3),znucl(ntypat),spinaxis(3)
2056 : : real(dp),intent(in) :: ratsph(ntypat),rhomag(2,nspden)
2057 : : real(dp),intent(in),optional :: ziontypat(ntypat)
2058 : :
2059 : : !Local variables ------------------------------
2060 : : !scalars
2061 : : integer :: iatom,icplex
2062 : : real(dp) :: mag_coll , mag_x, mag_y, mag_z ! EB
2063 : : real(dp) :: mag_coll_im, mag_x_im, mag_y_im, mag_z_im ! SPr
2064 : : real(dp) :: rho_tot, rho_tot_im
2065 : : real(dp) :: sum_mag, sum_mag_x,sum_mag_y,sum_mag_z,sum_rho_up,sum_rho_dn,sum_rho_tot ! EB
2066 : : real(dp) :: sum_mag_im, sum_mag_x_im,sum_mag_y_im,sum_mag_z_im,sum_rho_up_im,sum_rho_dn_im,sum_rho_tot_im ! SR
2067 : : real(dp) :: mag_r, mag_theta, mag_phi, vec(3)
2068 : : real(dp) :: sum_mag_r, sum_mag_theta, sum_mag_phi
2069 : : real(dp) :: exact_mag_r, exact_mag_theta, exact_mag_phi
2070 : : real(dp) :: exact_mag_r_im, exact_mag_theta_im, exact_mag_phi_im
2071 : : real(dp) :: alpha, beta, Rspin(3,3), Rspin_t(3,3), mag_cart(3), mag_spin(3), mag_tot_cart(3), mag_tot_cart_im(3)
2072 : : character(len=500) :: msg,msg1
2073 : : character(len=500) :: msg_cplex
2074 : : type(atomdata_t) :: atom
2075 : :
2076 : : ! *************************************************************************
2077 : : !write(ab_out,*)' prtdenmagsph : enter, rhomag(1,2)=',rhomag(1,2)
2078 : :
2079 [ + + ]: 9823 : if(nspden==2)then
2080 : 1326 : rho_tot=rhomag(1,1) ; mag_coll=rhomag(1,2)
2081 [ + + ]: 1326 : if(cplex==2)then
2082 : 12 : rho_tot_im=rhomag(2,1) ; mag_coll_im=rhomag(2,2)
2083 : : endif
2084 [ + + ]: 8497 : else if(nspden==4)then
2085 : 2146 : rho_tot=rhomag(1,1) ; mag_x=rhomag(1,2) ; mag_y=rhomag(1,3) ; mag_z=rhomag(1,4)
2086 [ + + ]: 2146 : if(cplex==2)then
2087 : 11 : rho_tot_im=rhomag(2,1) ; mag_x_im=rhomag(2,2) ; mag_y_im=rhomag(2,3) ; mag_z_im=rhomag(2,4)
2088 : : endif
2089 : : endif
2090 : :
2091 [ + - ]: 9823 : if(option/=0)then
2092 : :
2093 : : !Printing
2094 : 9823 : sum_mag=zero
2095 : 9823 : sum_mag_x=zero
2096 : 9823 : sum_mag_y=zero
2097 : 9823 : sum_mag_z=zero
2098 : 9823 : sum_mag_x_im=zero
2099 : 9823 : sum_mag_y_im=zero
2100 : 9823 : sum_mag_z_im=zero
2101 : 9823 : sum_rho_up=zero
2102 : 9823 : sum_rho_dn=zero
2103 : 9823 : sum_rho_tot=zero
2104 : :
2105 : 9823 : sum_mag_im=zero
2106 : 9823 : sum_mag_x_im=zero
2107 : 9823 : sum_mag_y_im=zero
2108 : 9823 : sum_mag_z_im=zero
2109 : 9823 : sum_rho_up_im=zero
2110 : 9823 : sum_rho_dn_im=zero
2111 : 9823 : sum_rho_tot_im=zero
2112 : :
2113 : : !Rotation matrices identity by default
2114 : 9823 : Rspin(:,:)=zero ; Rspin_t(:,:)=zero
2115 : 9823 : Rspin(1,1)=one ; Rspin(2,2)=one ; Rspin(3,3)=one
2116 : 9823 : Rspin_t(:,:)=Rspin(:,:)
2117 : :
2118 : : !Print spinaxis info only if axis not aligned with z
2119 [ + + - + ]: 9823 : if (abs(spinaxis(1))>tol8 .or. abs(spinaxis(2))>tol8) then
2120 : 1 : call geteuler(spinaxis, alpha, beta)
2121 : 1 : msg=' Spinaxis rotation information:'
2122 : 1 : write(msg, '(3a)' ) trim(msg),ch10,' ------------------------------'; call wrtout(units,msg)
2123 : 1 : write(msg, '(a,f12.6)') ' Alpha rotation angle around z-axis (degrees):', alpha * 180.0_dp / pi; call wrtout(units,msg)
2124 : 1 : write(msg, '(a,f12.6)') ' Beta rotation angle around y-axis (degrees): ', beta * 180.0_dp / pi; call wrtout(units,msg)
2125 : 1 : write(msg, '(a)') ' ---------------------------------------------------------'; call wrtout(units,msg)
2126 : 1 : call cart2spinaxis(alpha, beta, Rspin)
2127 : 1 : write(msg, '(a)') ' Rotation matrix from cartesian coordinate to spinaxis coordinate'; call wrtout(units,msg)
2128 : 1 : write(msg, '(3f14.8)') Rspin(1,1), Rspin(1,2), Rspin(1,3); call wrtout(units,msg)
2129 : 1 : write(msg, '(3f14.8)') Rspin(2,1), Rspin(2,2), Rspin(2,3); call wrtout(units,msg)
2130 : 1 : write(msg, '(3f14.8)') Rspin(3,1), Rspin(3,2), Rspin(3,3); call wrtout(units,msg)
2131 : :
2132 [ + + + + ]: 13 : Rspin_t = transpose(Rspin)
2133 : 1 : write(msg, '(a)') ' Rotation matrix from spinaxis coordinate to cartesian coordinate'; call wrtout(units,msg)
2134 : 1 : write(msg, '(3f14.8)') Rspin_t(1,1), Rspin_t(1,2), Rspin_t(1,3); call wrtout(units,msg)
2135 : 1 : write(msg, '(3f14.8)') Rspin_t(2,1), Rspin_t(2,2), Rspin_t(2,3); call wrtout(units,msg)
2136 : 1 : write(msg, '(3f14.8)') Rspin_t(3,1), Rspin_t(3,2), Rspin_t(3,3); call wrtout(units,msg)
2137 : 1 : write(msg, '(a)') ' ----------------------------------------------------------------'
2138 : 1 : call wrtout(units,msg)
2139 : : end if
2140 : :
2141 : 9823 : mag_tot_cart=zero
2142 : 9823 : mag_tot_cart_im=zero
2143 [ + + ]: 9823 : if (nspden==4) then
2144 [ + + ]: 8584 : mag_spin = [mag_x,mag_y,mag_z]
2145 [ + + + + ]: 27898 : mag_tot_cart = matmul(Rspin_t,mag_spin)
2146 [ + + ]: 2146 : if (cplex==2) then
2147 [ + + ]: 44 : mag_spin = [mag_x_im,mag_y_im,mag_z_im]
2148 [ + + + + ]: 143 : mag_tot_cart_im = matmul(Rspin_t,mag_spin)
2149 : : end if
2150 : : end if
2151 : :
2152 [ + - - - ]: 9823 : if(option==1 .or. option==11 .or. option==21) then
2153 : :
2154 [ + + ]: 9823 : if(nspden==1) then
2155 [ + + ]: 6351 : if(option== 1)msg1=' Integrated electronic density in atomic spheres:'
2156 [ + + ]: 6351 : if(option==11)msg1=ch10//' Integrated potential residual in atomic spheres:'
2157 [ + + ]: 6351 : if(option==21)msg1=ch10//' Gradient with respect to target (=torque) :'
2158 : 6351 : write(msg, '(3a)' ) trim(msg1),ch10,' ------------------------------------------------'
2159 : 6351 : call wrtout(units,msg)
2160 [ + + ]: 6351 : if(ratsm>tol8)then
2161 : 221 : write(msg, '(a,f8.4,a)' ) ' Radius=ratsph(iatom), smearing ratsm=',ratsm,'.'
2162 : 221 : call wrtout(units,msg)
2163 : : endif
2164 [ + + ]: 6351 : if(option== 1)then
2165 : 6241 : msg=' Atom Sphere_radius Integrated_density'
2166 [ + + ]: 6241 : if(present(ziontypat)) write(msg,'(a,a)')trim(msg),' Atomic charge'
2167 : : endif
2168 [ + + ]: 6351 : if(option==11)msg=' Atom Sphere_radius Integrated_potresid'
2169 [ + + ]: 6351 : if(option==21)msg=' Atom Sphere_radius Torque'
2170 : 6351 : call wrtout(units,msg)
2171 [ + + ]: 25576 : do iatom=1,natom
2172 : 19225 : call atomdata_from_znucl(atom, znucl(typat(iatom)))
2173 : 19225 : write(msg, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),intgden(1,1,iatom)
2174 : : !write(msg, '(i5,a3,f15.5,f20.8)' ) iatom,atom%symbol,ratsph(typat(iatom)),intgden(1,1,iatom)
2175 [ + + ]: 19225 : if(option==21)then
2176 : : ! There is a change of sign to get the gradient wrt chrgat.
2177 : 10 : write(msg, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),-intgden(1,1,iatom)
2178 : : !write(msg, '(i5,a3,f15.5,f20.8)' ) iatom,atom%symbol,ratsph(typat(iatom)),-intgden(1,1,iatom)
2179 : : endif
2180 : : !If option=1, print atomic charge
2181 [ + + + + ]: 19225 : if(option==1 .and. present(ziontypat))then
2182 : 10 : write(msg, '(a,f20.8)' ) trim(msg),ziontypat(typat(iatom))-intgden(1,1,iatom)
2183 : : endif
2184 : 44801 : call wrtout(units,msg)
2185 : : end do
2186 : : endif
2187 : :
2188 [ + + ]: 9823 : if(nspden==2 .or. nspden==4) then
2189 : :
2190 [ + + ]: 3472 : if(option== 1)msg1=' Integrated electronic and magnetization densities in atomic spheres:'
2191 [ + + ]: 3472 : if(option==11)msg1=ch10//' Integrated potential residual in atomic spheres (scalar + magnetic field):'
2192 [ + + ]: 3472 : if(option==21)msg1=ch10//' Gradient with respect to target (=torque) (scalar + magnetic field):'
2193 : 3472 : write(msg, '(3a)' ) trim(msg1),ch10,' ------------------------------------------------------------------------------------------------------------'
2194 : 3472 : call wrtout(units,msg)
2195 : :
2196 [ + + + + ]: 3472 : if(option== 1 .and. nspden==2) msg1='. Diff(up-dn)=approximate z local magnetic moment.'
2197 [ + + + + ]: 3472 : if(option== 1 .and. nspden==4) msg1='. mag(i)=approximate local magnetic moment.'
2198 [ + + ]: 3472 : if(option==11 .or. option==21) msg1='.'
2199 : 3472 : write(msg, '(a,f8.4,a)' ) ' Radius=ratsph(iatom), smearing ratsm=',ratsm,trim(msg1)
2200 : 3472 : call wrtout(units,msg)
2201 : :
2202 [ + + ]: 3472 : if(option==1)then
2203 [ + + ]: 2234 : if(nspden==2) msg=' Atom Radius up_density dn_density Total(up+dn) Diff(up-dn)'
2204 [ + + ]: 2234 : if(nspden==4) msg=' Atom Radius Total density |mag| mag(x) mag(y) mag(z) mag(theta) mag(phi) '
2205 [ + + + - ]: 2234 : if(present(ziontypat))msg=trim(msg)//' Atomic charge'
2206 [ + + ]: 1238 : else if(option==11)then
2207 [ + + ]: 1189 : if(nspden==2) msg=' Atom Radius Potential B(z) up pot down pot'
2208 [ + + ]: 1189 : if(nspden==4) msg=' Atom Radius Potential |B| B(x) B(y) B(z) B(theta) B(phi) '
2209 [ + - ]: 49 : else if(option==21)then
2210 [ + + ]: 49 : if(nspden==2) msg=' Atom Radius grchrg T(z) up torque down torque'
2211 [ + + ]: 49 : if(nspden==4) msg=' Atom Radius Torque |T| T(x) T(y) T(z) T(theta) T(phi) '
2212 : : endif
2213 [ + + ]: 3472 : if (cplex==2) then
2214 : 23 : msg_cplex=msg
2215 : : else
2216 : 3449 : call wrtout(units,msg)
2217 : : endif
2218 : :
2219 : : endif
2220 : :
2221 [ + + ]: 9823 : if(nspden==2)then
2222 [ + + ]: 2664 : do icplex=1,cplex
2223 : 1338 : sum_mag=zero
2224 : 1338 : sum_rho_up=zero
2225 : 1338 : sum_rho_dn=zero
2226 : 1338 : sum_rho_tot=zero
2227 [ + + ]: 1338 : if (cplex==2) then
2228 [ + + ]: 24 : if (icplex==1) then
2229 : 12 : write(msg, '(a)') ' Real part of magnetization:'
2230 : 12 : call wrtout(units,msg)
2231 : 12 : call wrtout(units,msg_cplex)
2232 [ + - ]: 12 : elseif (icplex==2) then
2233 : 12 : write(msg, '(a)') '----------------------------'
2234 : 12 : call wrtout(units,msg)
2235 : 12 : write(msg, '(a)') ' Imaginary part of magnetization:'
2236 : 12 : call wrtout(units,msg)
2237 : 12 : call wrtout(units,msg_cplex)
2238 : : endif
2239 : : endif
2240 [ + + ]: 3854 : do iatom=1,natom
2241 : 2516 : call atomdata_from_znucl(atom, znucl(typat(iatom)))
2242 [ + + ]: 2516 : if(option/=21)then
2243 : 2492 : write(msg,'(i5,a3,f10.5,2f13.6,a,f12.6,a,f12.6)' ) iatom,atom%symbol,ratsph(typat(iatom)),intgden(icplex,1,iatom),intgden(icplex,2,iatom)
2244 : : else
2245 : 24 : write(msg,'(i5,a3,f10.5,2f13.6,a,f12.6,a,f12.6)' ) iatom,atom%symbol,ratsph(typat(iatom)),-intgden(icplex,1,iatom),intgden(icplex,2,iatom)
2246 : : endif
2247 : 2516 : write(msg,'(a,a,f12.6,a,f12.6)')trim(msg),' ',(intgden(icplex,1,iatom)+intgden(icplex,2,iatom)),' ',(intgden(icplex,1,iatom)-intgden(icplex,2,iatom))
2248 [ + + + + ]: 2516 : if(option==1 .and. present(ziontypat)) then
2249 : 24 : write(msg, '(a,f14.6)') trim(msg),ziontypat(typat(iatom))-(intgden(icplex,1,iatom)+intgden(icplex,2,iatom))
2250 : : end if
2251 : 2516 : call wrtout(units,msg)
2252 : : ! Compute the sum of the magnetization
2253 : 2516 : sum_mag=sum_mag+intgden(icplex,1,iatom)-intgden(icplex,2,iatom)
2254 : 2516 : sum_rho_up=sum_rho_up+intgden(icplex,1,iatom)
2255 : 2516 : sum_rho_dn=sum_rho_dn+intgden(icplex,2,iatom)
2256 : 6370 : sum_rho_tot=sum_rho_tot+intgden(icplex,1,iatom)+intgden(icplex,2,iatom)
2257 : : end do
2258 : 1338 : write(msg, '(a)') ' ------------------------------------------------------------------------------------------------------------'
2259 : 1338 : call wrtout(units,msg)
2260 : 1338 : write(msg, '(a,2f13.6,a,f12.6,a,f12.6)') ' Sum: ', sum_rho_up,sum_rho_dn,' ',sum_rho_tot,' ',sum_mag
2261 : 1338 : call wrtout(units,msg)
2262 : :
2263 [ + + ]: 1338 : if(option==1)then
2264 : 1071 : write(msg, '(a,f14.6)') ' Total magnetization (from the atomic spheres): ', sum_mag
2265 : 1071 : call wrtout(units,msg)
2266 [ + + ]: 1071 : if (icplex==1) then
2267 : 1059 : write(msg, '(a,f14.6)') ' Total magnetization (exact up - dn): ', mag_coll
2268 [ + - ]: 12 : elseif (icplex==2) then
2269 : 12 : write(msg, '(a,f14.6)') ' Total magnetization (exact up - dn): ', mag_coll_im
2270 : : endif
2271 : 1071 : call wrtout(units,msg)
2272 : : endif
2273 : 1338 : write(msg, '(a)') ' '
2274 : 2664 : call wrtout(units,msg)
2275 : : enddo
2276 : :
2277 [ + + ]: 8497 : elseif(nspden==4) then
2278 : :
2279 [ + + ]: 4303 : do icplex=1,cplex
2280 : 2157 : sum_mag_x=zero
2281 : 2157 : sum_mag_y=zero
2282 : 2157 : sum_mag_z=zero
2283 [ + + ]: 2157 : if (cplex==2) then
2284 [ + + ]: 22 : if (icplex==1) then
2285 : 11 : write(msg, '(a)') ' Real part of magnetization:'
2286 : 11 : call wrtout(units,msg)
2287 : 11 : call wrtout(units,msg_cplex)
2288 [ + - ]: 11 : elseif (icplex==2) then
2289 : 11 : write(msg, '(a)') '----------------------------'
2290 : 11 : call wrtout(units,msg)
2291 : 11 : write(msg, '(a)') ' Imaginary part of magnetization:'
2292 : 11 : call wrtout(units,msg)
2293 : 11 : call wrtout(units,msg_cplex)
2294 : : endif
2295 : : endif
2296 [ + + ]: 7633 : do iatom=1,natom
2297 [ + + + + ]: 71188 : mag_cart(:)=matmul(Rspin_t,intgden(icplex,2:4,iatom))
2298 : 5476 : call vcart2ylm(mag_cart, mag_r, mag_theta, mag_phi)
2299 : 5476 : call atomdata_from_znucl(atom, znucl(typat(iatom)))
2300 [ + + ]: 5476 : if(option/=21)then
2301 : 5384 : write(msg, '(i5,a3,f10.5,f16.6,a,6f12.6)' ) iatom,atom%symbol,ratsph(typat(iatom)),intgden(icplex,1,iatom),' ',mag_r,mag_cart(1),mag_cart(2),mag_cart(3),mag_theta,mag_phi
2302 : : else
2303 : 92 : write(msg, '(i5,a3,f10.5,f16.6,a,6f12.6)' ) iatom,atom%symbol,ratsph(typat(iatom)),-intgden(icplex,1,iatom),' ',mag_r,mag_cart(1),mag_cart(2),mag_cart(3),mag_theta,mag_phi
2304 : : endif
2305 [ + + + + ]: 5476 : if(option==1 .and. present(ziontypat))&
2306 : 92 : & write(msg, '(a,f14.6)') trim(msg),ziontypat(typat(iatom))-intgden(icplex,1,iatom)
2307 : 5476 : call wrtout(units,msg)
2308 : : ! Compute the sum of the magnetization in x, y and z directions
2309 : 5476 : sum_mag_x=sum_mag_x+mag_cart(1)
2310 : 5476 : sum_mag_y=sum_mag_y+mag_cart(2)
2311 : 13109 : sum_mag_z=sum_mag_z+mag_cart(3)
2312 : : enddo
2313 [ + + ]: 8628 : vec = (/sum_mag_x,sum_mag_y,sum_mag_z /)
2314 : 2157 : call vcart2ylm(vec, sum_mag_r, sum_mag_theta, sum_mag_phi)
2315 : 2157 : write(msg, '(a)') ' ------------------------------------------------------------------------------------------------------------'
2316 : 2157 : call wrtout(units,msg)
2317 : :
2318 [ + + ]: 6460 : if(option==1)then
2319 [ + + + + : 4710 : if (any(qgbt /= zero)) then
+ + ]
2320 : 14 : write(msg, '(a,f12.6,f12.6,f12.6)') ' Cell sum of sphere magnetization'
2321 : 14 : call wrtout(units,msg)
2322 : : else
2323 : 1172 : write(msg, '(a)') ' Total magnetization (spheres) '
2324 : 1172 : call wrtout(units,msg)
2325 : : endif
2326 : 1186 : write(msg, '(a,3f12.6)') ' (cart.coord.) ', sum_mag_x,sum_mag_y,sum_mag_z
2327 : 1186 : call wrtout(units,msg)
2328 : 1186 : write(msg, '(a,f12.6,a,2f12.6)') ' (sph.coord.) ', sum_mag_r," ",sum_mag_theta,sum_mag_phi
2329 : 1186 : call wrtout(units,msg)
2330 : :
2331 [ + + + + : 4710 : if (any(qgbt /= zero)) then
+ + ]
2332 : 14 : write(msg, '(a,f10.6,f12.6,f12.6)') ' Integral of periodic magnetization'
2333 : 14 : call wrtout(units,msg)
2334 : : else
2335 : 1172 : write(msg, '(a,f12.6,f12.6,4f12.6)') ' Total magnetization (exact) '
2336 : 1172 : call wrtout(units,msg)
2337 : : endif
2338 [ + + ]: 1186 : if (icplex==1) then
2339 : 1175 : vec = mag_tot_cart
2340 : 1175 : call vcart2ylm(vec, exact_mag_r, exact_mag_theta, exact_mag_phi)
2341 : 1175 : write(msg, '(a,3f12.6)') ' (cart.coord.) ', mag_tot_cart(1),mag_tot_cart(2),mag_tot_cart(3)
2342 : 1175 : call wrtout(units,msg)
2343 : 1175 : write(msg, '(a,f12.6,a,2f12.6)') ' (sph.coord.) ', exact_mag_r," ",exact_mag_theta,exact_mag_phi
2344 : 1175 : call wrtout(units,msg)
2345 [ + - ]: 11 : elseif (icplex==2) then
2346 : 11 : vec = mag_tot_cart_im
2347 : 11 : call vcart2ylm(vec, exact_mag_r_im, exact_mag_theta_im, exact_mag_phi_im)
2348 : 11 : write(msg, '(a,3f12.6)') ' (cart.coord.) ', mag_tot_cart_im(1),mag_tot_cart_im(2),mag_tot_cart_im(3)
2349 : 11 : call wrtout(units,msg)
2350 : 11 : write(msg,'(a,f12.6,a,2f12.6)')' (sph.coord.) ', exact_mag_r_im," ", exact_mag_theta_im,exact_mag_phi_im
2351 : 11 : call wrtout(units,msg)
2352 : : endif
2353 : : endif
2354 : : enddo
2355 : : end if
2356 : :
2357 : : elseif(option==-1) then
2358 : :
2359 : 0 : write(msg, '(2a)') ch10,' ---------------------------------------------------------------------------------------------------------------'
2360 : 0 : call wrtout(units,msg)
2361 : :
2362 [ # # ]: 0 : if(nspden==1) then
2363 : : write(msg, '(4a)' ) &
2364 : 0 : & ' Fermi level charge density n_f:',ch10,&
2365 : 0 : & ' ------------------------------------------------------------------------',ch10
2366 : : else
2367 : : write(msg, '(4a)' ) &
2368 : 0 : & ' Fermi level charge density n_f and magnetization m_f:',ch10,&
2369 : 0 : & ' ------------------------------------------------------------------------',ch10
2370 : : end if
2371 : 0 : call wrtout(units,msg)
2372 : :
2373 [ # # ]: 0 : if(cplex==1) then
2374 : 0 : write(msg, '(a,f13.8)') ' n_f = ',rho_tot
2375 : : else
2376 : 0 : write(msg, '(a,f13.8,a,f13.8)') ' Re[n_f]= ', rho_tot," Im[n_f]= ",rho_tot_im
2377 : : end if
2378 : 0 : call wrtout(units,msg)
2379 [ # # ]: 0 : if(nspden==2) then
2380 [ # # ]: 0 : if(cplex==1) then
2381 : 0 : write(msg, '(a,f13.8)') ' m_f = ', mag_coll
2382 : : else
2383 : 0 : write(msg, '(a,f13.8,a,f13.8)') ' Re[m_f]= ', mag_coll," Im[m_f]= ",mag_coll_im
2384 : : end if
2385 : 0 : call wrtout(units,msg)
2386 [ # # ]: 0 : elseif (nspden==4) then
2387 : 0 : write(msg, '(a,f13.8)') ' mx_f = ',mag_tot_cart(1)
2388 : 0 : call wrtout(units,msg)
2389 : 0 : write(msg, '(a,f13.8)') ' my_f = ',mag_tot_cart(2)
2390 : 0 : call wrtout(units,msg)
2391 : 0 : write(msg, '(a,f13.8)') ' mz_f = ',mag_tot_cart(3)
2392 : 0 : call wrtout(units,msg)
2393 : : end if
2394 : :
2395 : 0 : write(msg, '(3a)') ch10,' ------------------------------------------------------------------------',ch10
2396 : 0 : call wrtout(units,msg)
2397 : :
2398 : :
2399 : : else if (option==2 .or. option==3 .or. option==4) then
2400 : : ! Used in the DFPT case, option=idir+1
2401 : :
2402 [ # # ]: 0 : if(abs(rho_tot)<1.0d-10) then
2403 : 0 : rho_tot=0
2404 : : end if
2405 : :
2406 : 0 : write(msg, '(2a)') ch10,' ---------------------------------------------------------------------------'
2407 : 0 : call wrtout(units,msg)
2408 : :
2409 : : ! if(nspden==1) then
2410 : : ! write(msg, '(4a)' ) &
2411 : : !& ' Integral of the first order density n^(1):',ch10,&
2412 : : !& ' ------------------------------------------------------------------------',ch10
2413 : : ! else
2414 : : ! write(msg, '(4a)' ) &
2415 : : !& ' Integrals of the first order density n^(1) and magnetization m^(1):',ch10,&
2416 : : !& ' ------------------------------------------------------------------------',ch10
2417 : : ! end if
2418 : : ! call wrtout(units,msg)
2419 : : !
2420 : : ! if(cplex==1) then
2421 : : ! write(msg, '(a,e16.8)') ' n^(1) = ', rho_tot
2422 : : ! else
2423 : : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[n^(1)] = ', rho_tot," Im[n^(1)] = ",rho_tot_im
2424 : : ! end if
2425 : : ! call wrtout(units,msg)
2426 : : !
2427 : : ! if(nspden==2) then
2428 : : !
2429 : : ! if(cplex==1) then
2430 : : ! write(msg, '(a,e16.8)') ' m^(1) = ', mag_coll
2431 : : ! else
2432 : : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[m^(1)] = ', mag_coll," Im[m^(1)] = ",mag_coll_im
2433 : : ! end if
2434 : : ! call wrtout(units,msg)
2435 : : !
2436 : : ! elseif (nspden==4) then
2437 : : ! if(cplex==1) then
2438 : : ! write(msg, '(a,e16.8)') ' mx^(1) = ', mag_x
2439 : : ! call wrtout(units,msg)
2440 : : ! write(msg, '(a,e16.8)') ' my^(1) = ', mag_y
2441 : : ! call wrtout(units,msg)
2442 : : ! write(msg, '(a,e16.8)') ' mz^(1) = ', mag_z
2443 : : ! call wrtout(units,msg)
2444 : : ! else
2445 : : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[mx^(1)]= ', mag_x, " Im[mx^(1)]= ", mag_x_im
2446 : : ! call wrtout(units,msg)
2447 : : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[my^(1)]= ', mag_y, " Im[my^(1)]= ", mag_y_im
2448 : : ! call wrtout(units,msg)
2449 : : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[mz^(1)]= ', mag_z, " Im[mz^(1)]= ", mag_z_im
2450 : : ! call wrtout(units,msg)
2451 : : ! end if
2452 : : ! end if
2453 : : !
2454 : : ! write(msg, '(3a)') ch10,' ------------------------------------------------------------------------',ch10
2455 : : ! call wrtout(units,msg)
2456 : :
2457 : : end if
2458 : :
2459 : : end if ! option/=0
2460 : :
2461 : 9823 : end subroutine prtdenmagsph
2462 : : !!***
2463 : :
2464 : : !!****f* m_dens/radsmear
2465 : : !! NAME
2466 : : !! radsmear
2467 : : !!
2468 : : !! FUNCTION
2469 : : !! As a function of the argument xarg (a positive number), return a function fsm that is zero
2470 : : !! beyond some cut-off value xcut, one for xarg smaller than xcut-xsmear,
2471 : : !! and interpolates smoothly between one and zero in the region from xcut-xsmear to xcut.
2472 : : !! Also returns the derivative of this function, called dfsm.
2473 : : !! The function fsm is twice differentiable at xcut (the first derivative is continuous, not the second),
2474 : : !! and three times differentiable at xcut-xsmear (the second derivative is continuous, not the third).
2475 : : !!
2476 : : !!
2477 : : !! INPUTS
2478 : : !! xarg=argument of the function (should be positive)
2479 : : !! xcut=largest value for which the function is non-zero
2480 : : !! xsmear=defined the smearing region, between xcut-xsmear and xcut
2481 : : !!
2482 : : !! OUTPUT
2483 : : !! fsm=value of the function
2484 : : !! dfsm=derivative of the function with respect to xarg (zero, except in the smearing region).
2485 : : !!
2486 : : !! SOURCE
2487 : :
2488 : 97135536 : subroutine radsmear(dfsm,fsm,xarg,xcut,xsmear)
2489 : :
2490 : : !Arguments ------------------------------------
2491 : : !scalars
2492 : : real(dp), intent(out) :: dfsm,fsm
2493 : : real(dp), intent(in) :: xarg, xcut, xsmear
2494 : :
2495 : : !Local variables ------------------------------
2496 : : !scalars
2497 : : real(dp) :: xsmearinv,xx
2498 : : !******************************************************************
2499 : :
2500 : 97135536 : fsm = zero
2501 : 97135536 : dfsm=zero
2502 [ + + ]: 97135536 : if (xarg < xcut - xsmear - tol12) then
2503 : 95944818 : fsm = one
2504 [ + - ]: 1190718 : else if (xarg < xcut - tol12) then
2505 : 1190718 : xsmearinv=one/xsmear
2506 : 1190718 : xx = (xcut - xarg) * xsmearinv
2507 : 1190718 : fsm = xx**2*(3+xx*(1+xx*(-6+3*xx)))
2508 : 1190718 : dfsm = -(xx*(6+xx*(3+xx*(-24+15*xx))))*xsmearinv
2509 : : end if
2510 : :
2511 : 97135536 : end subroutine radsmear
2512 : : !!***
2513 : :
2514 : : !!****f* ABINIT/printmagvtk
2515 : : !! NAME
2516 : : !! printmagvtk
2517 : : !!
2518 : : !! FUNCTION
2519 : : !! Auxiliary routine for printing out magnetization density in VTK format.
2520 : : !! Output file name is DEN.vtk
2521 : : !!
2522 : : !! INPUTS
2523 : : !! mpi_enreg = information about adopted parallelization strategy
2524 : : !! nspden = number of components of density matrix (possible values re 1,2, or 4)
2525 : : !! nspden: 1 -> rho
2526 : : !! nspden: 2 -> rho_up,rho_dwn
2527 : : !! nspden: 4 -> rho,mx,my,mz
2528 : : !! nfft = number of fft points per FFT processor
2529 : : !! ngfft = full information about FFT mesh
2530 : : !! rhor = density array stored in the memory of current FFT processor
2531 : : !! rprimd = array of lattice vectors
2532 : : !!
2533 : : !! OUTPUT
2534 : : !!
2535 : : !! NOTES
2536 : : !! At the moment this routine is mainly used for development and debugging
2537 : : !! of gs and dfpt calculations with non-collinear spins. If needed, can be used
2538 : : !! to print final density in vtk format.
2539 : : !! IMPORTANT: implementation is thoroughly checked only for npspinor = 1,
2540 : : !! for other case might need to change the part gathering
2541 : : !! the FFT mesh info
2542 : : !!
2543 : : !! SOURCE
2544 : :
2545 : : subroutine printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,fname)
2546 : :
2547 : : !Arguments ------------------------------------
2548 : : !scalars
2549 : : type(MPI_type),intent(in) :: mpi_enreg
2550 : : integer,intent(in) :: nfft,nspden,cplex
2551 : : !arrays
2552 : : integer,intent(in) :: ngfft(18)
2553 : : real(dp),intent(in) :: rhor(cplex*nfft,nspden),rprimd(3,3)
2554 : : character(len=*),intent(in) :: fname
2555 : :
2556 : : !Local variables-------------------------------
2557 : : !scalars
2558 : : integer :: denvtk,denxyz,denxyz_im,nfields
2559 : : integer :: nx,ny,nz,nfft_tot
2560 : : integer :: ii,jj,kk,ind,ispden
2561 : : integer :: mpi_comm,mpi_head,mpi_rank,ierr
2562 : : real(dp) :: rx,ry,rz
2563 : : integer :: nproc_fft,ir
2564 : : character(len=500) :: msg
2565 : : character(len=10) :: outformat
2566 : : character(len=fnlen) :: fname_vtk, fname_xyz, fname_xyz_re, fname_xyz_im
2567 : : !arrays
2568 : : real(dp),allocatable :: rhorfull(:,:)
2569 : : ! *************************************************************************
2570 : :
2571 : : fname_vtk=adjustl(adjustr(fname)//".vtk")
2572 : : fname_xyz=adjustl(adjustr(fname)//".xyz")
2573 : : fname_xyz_re=adjustl(adjustr(fname)//"_re.xyz")
2574 : : fname_xyz_im=adjustl(adjustr(fname)//"_im.xyz")
2575 : : !write(std_out,*) ' Writing out .vtk file: ',fname_vtk
2576 : : !write(std_out,*) ' Writing out .xyz file: ',fname_xyz
2577 : :
2578 : : !if 1 or two component density then write out either 1 or 2 scalar density fields
2579 : : !if 4, then write one scalar field (density) and one vector field (magnetization density)
2580 : : if(nspden/=4)then
2581 : : nfields=nspden
2582 : : else
2583 : : nfields=2
2584 : : end if
2585 : :
2586 : : nfields=nfields*cplex
2587 : :
2588 : : ! FFT mesh specifications: full grid
2589 : : nx=ngfft(1) ! number of points along 1st lattice vector
2590 : : ny=ngfft(2) ! number of points along 2nd lattice vector
2591 : : nz=ngfft(3) ! number of points along 3rd lattice vector
2592 : : nfft_tot=nx*ny*nz ! total number of fft mesh points (can be different from nfft in case of distributed memory of nproc_fft processors)
2593 : :
2594 : :
2595 : : ! Gather information about memory distribution
2596 : : mpi_head=0
2597 : : mpi_comm = mpi_enreg%comm_fft
2598 : : mpi_rank = xmpi_comm_rank(mpi_comm)
2599 : : nproc_fft=ngfft(10)
2600 : :
2601 : : ! Create array to host full FFT mesh
2602 : : if(mpi_rank==mpi_head)then
2603 : : ABI_MALLOC(rhorfull,(cplex*nfft_tot,nspden))
2604 : : end if
2605 : :
2606 : : ! Fill in the full mesh
2607 : : if(nproc_fft==1)then
2608 : : rhorfull=rhor
2609 : : else
2610 : : do ir=1,nspden
2611 : : call xmpi_gather(rhor(:,ir),cplex*nfft,rhorfull(:,ir),cplex*nfft,mpi_head,mpi_comm,ierr)
2612 : : end do
2613 : : end if
2614 : :
2615 : : if(mpi_rank==mpi_head)then
2616 : :
2617 : : ! Open the output vtk file
2618 : : if (open_file(fname_vtk,msg,newunit=denvtk,status='replace',form='formatted') /=0) then
2619 : : ABI_WARNING(msg)
2620 : : RETURN
2621 : : end if
2622 : :
2623 : : if(cplex==1) then
2624 : : if (open_file(fname_xyz,msg,newunit=denxyz,status='replace',form='formatted') /=0) then
2625 : : ABI_WARNING(msg)
2626 : : RETURN
2627 : : end if
2628 : : else if (cplex==2) then
2629 : : if (open_file(fname_xyz_re,msg,newunit=denxyz,status='replace',form='formatted') /=0) then
2630 : : ABI_WARNING(msg)
2631 : : RETURN
2632 : : end if
2633 : : if (open_file(fname_xyz_im,msg,newunit=denxyz_im,status='replace',form='formatted') /=0) then
2634 : : ABI_WARNING(msg)
2635 : : RETURN
2636 : : end if
2637 : : end if
2638 : :
2639 : : ! Write the header of the output vtk file
2640 : : write(denvtk,"(a)") '# vtk DataFile Version 2.0'
2641 : : write(denvtk,"(a)") 'Electron density components'
2642 : : write(denvtk,"(a)") 'ASCII'
2643 : : write(denvtk,"(a)") 'DATASET STRUCTURED_GRID'
2644 : : write(denvtk,"(a,3i6)") 'DIMENSIONS ', nx,ny,nz
2645 : : write(denvtk,"(a,i18,a)") 'POINTS ',nfft_tot,' double'
2646 : :
2647 : : if (nspden==1) then
2648 : : outformat="(4e16.8)"
2649 : : else if (nspden==2) then
2650 : : outformat="(5e16.8)"
2651 : : else
2652 : : outformat="(7e16.8)"
2653 : : end if
2654 : :
2655 : : ! Write out information about grid points
2656 : : do kk=0,nz-1
2657 : : do jj=0,ny-1
2658 : : do ii=0,nx-1
2659 : :
2660 : : rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3)
2661 : : ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3)
2662 : : rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3)
2663 : : write(denvtk,'(3f16.8)') rx,ry,rz !coordinates of the grid point
2664 : : ind=1+ii+nx*(jj+ny*kk)
2665 : : if (cplex==1) then
2666 : : write(denxyz,outformat) rx,ry,rz,(rhorfull(ind,ispden),ispden=1,nspden)
2667 : : else
2668 : : write(denxyz,outformat) rx,ry,rz,(rhorfull(2*ind-1,ispden),ispden=1,nspden)
2669 : : write(denxyz_im,outformat) rx,ry,rz,(rhorfull(2*ind ,ispden),ispden=1,nspden)
2670 : : end if
2671 : : end do
2672 : : end do
2673 : : end do
2674 : :
2675 : : if(cplex==1) then
2676 : : close(denxyz)
2677 : : else
2678 : : close(denxyz)
2679 : : close(denxyz_im)
2680 : : end if
2681 : :
2682 : : ! Write out information about field defined on the FFT mesh
2683 : : write(denvtk,"(a,i18)") 'POINT_DATA ',nfft_tot
2684 : : write(denvtk,"(a,i6)") 'FIELD Densities ',nfields
2685 : :
2686 : :
2687 : : ! Write out different fields depending on the number of density matrix components
2688 : : if(nspden==1)then
2689 : :
2690 : : !single component, so just write out the density
2691 : : if(cplex==1) then
2692 : : write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
2693 : : do kk=0,nz-1
2694 : : do jj=0,ny-1
2695 : : do ii=0,nx-1
2696 : : ind=1+ii+nx*(jj+ny*kk)
2697 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2698 : : end do
2699 : : end do
2700 : : end do
2701 : : else
2702 : : write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
2703 : : do kk=0,nz-1
2704 : : do jj=0,ny-1
2705 : : do ii=0,nx-1
2706 : : ind=2*(1+ii+nx*(jj+ny*kk))-1
2707 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2708 : : end do
2709 : : end do
2710 : : end do
2711 : : write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double'
2712 : : do kk=0,nz-1
2713 : : do jj=0,ny-1
2714 : : do ii=0,nx-1
2715 : : ind=2*(1+ii+nx*(jj+ny*kk))
2716 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2717 : : end do
2718 : : end do
2719 : : end do
2720 : : end if
2721 : :
2722 : : else if(nspden==2)then
2723 : :
2724 : : !two component, write the density for spin_up and spin_down channels
2725 : : if(cplex==1) then
2726 : :
2727 : : write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
2728 : : do kk=0,nz-1
2729 : : do jj=0,ny-1
2730 : : do ii=0,nx-1
2731 : : ind=1+ii+nx*(jj+ny*kk)
2732 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2733 : : end do
2734 : : end do
2735 : : end do
2736 : : write(denvtk,"(a,i18,a)") 'mag 1 ',nfft_tot,' double'
2737 : : do kk=0,nz-1
2738 : : do jj=0,ny-1
2739 : : do ii=0,nx-1
2740 : : ind=1+ii+nx*(jj+ny*kk)
2741 : : write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
2742 : : end do
2743 : : end do
2744 : : end do
2745 : :
2746 : : else
2747 : :
2748 : : write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
2749 : : do kk=0,nz-1
2750 : : do jj=0,ny-1
2751 : : do ii=0,nx-1
2752 : : ind=2*(1+ii+nx*(jj+ny*kk))-1
2753 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2754 : : end do
2755 : : end do
2756 : : end do
2757 : : write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double'
2758 : : do kk=0,nz-1
2759 : : do jj=0,ny-1
2760 : : do ii=0,nx-1
2761 : : ind=2*(1+ii+nx*(jj+ny*kk))
2762 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2763 : : end do
2764 : : end do
2765 : : end do
2766 : : write(denvtk,"(a,i18,a)") 'Re_mag 1 ',nfft_tot,' double'
2767 : : do kk=0,nz-1
2768 : : do jj=0,ny-1
2769 : : do ii=0,nx-1
2770 : : ind=2*(1+ii+nx*(jj+ny*kk))-1
2771 : : write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
2772 : : end do
2773 : : end do
2774 : : end do
2775 : : write(denvtk,"(a,i18,a)") 'Im_mag 1 ',nfft_tot,' double'
2776 : : do kk=0,nz-1
2777 : : do jj=0,ny-1
2778 : : do ii=0,nx-1
2779 : : ind=2*(1+ii+nx*(jj+ny*kk))
2780 : : write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
2781 : : end do
2782 : : end do
2783 : : end do
2784 : :
2785 : : end if
2786 : :
2787 : : else !here is the last option: nspden==4
2788 : :
2789 : : if(cplex==1) then
2790 : :
2791 : : !four component, write the density (scalar field) and magnetization density (vector field)
2792 : : write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
2793 : : do kk=0,nz-1
2794 : : do jj=0,ny-1
2795 : : do ii=0,nx-1
2796 : : ind=1+ii+nx*(jj+ny*kk)
2797 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2798 : : end do
2799 : : end do
2800 : : end do
2801 : : write(denvtk,"(a,i18,a)") 'mag 3 ',nfft_tot,' double'
2802 : : do kk=0,nz-1
2803 : : do jj=0,ny-1
2804 : : do ii=0,nx-1
2805 : : ind=1+ii+nx*(jj+ny*kk)
2806 : : write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
2807 : : end do
2808 : : end do
2809 : : end do
2810 : :
2811 : : else
2812 : :
2813 : : write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
2814 : : do kk=0,nz-1
2815 : : do jj=0,ny-1
2816 : : do ii=0,nx-1
2817 : : ind=2*(1+ii+nx*(jj+ny*kk))-1
2818 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2819 : : end do
2820 : : end do
2821 : : end do
2822 : : write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double'
2823 : : do kk=0,nz-1
2824 : : do jj=0,ny-1
2825 : : do ii=0,nx-1
2826 : : ind=2*(1+ii+nx*(jj+ny*kk))
2827 : : write(denvtk,'(f16.8)') rhorfull(ind,1)
2828 : : end do
2829 : : end do
2830 : : end do
2831 : : write(denvtk,"(a,i18,a)") 'Re_mag 3 ',nfft_tot,' double'
2832 : : do kk=0,nz-1
2833 : : do jj=0,ny-1
2834 : : do ii=0,nx-1
2835 : : ind=2*(1+ii+nx*(jj+ny*kk))-1
2836 : : write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
2837 : : end do
2838 : : end do
2839 : : end do
2840 : : write(denvtk,"(a,i18,a)") 'Im_mag 3 ',nfft_tot,' double'
2841 : : do kk=0,nz-1
2842 : : do jj=0,ny-1
2843 : : do ii=0,nx-1
2844 : : ind=2*(1+ii+nx*(jj+ny*kk))
2845 : : write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
2846 : : end do
2847 : : end do
2848 : : end do
2849 : :
2850 : : end if
2851 : :
2852 : : end if ! nspden options condition
2853 : :
2854 : : close (denvtk)
2855 : :
2856 : : !clean up the gathered FFT mesh
2857 : : ABI_FREE(rhorfull)
2858 : :
2859 : : end if
2860 : :
2861 : : end subroutine printmagvtk
2862 : : !!***
2863 : :
2864 : : !!****f* ABINIT/calmaxdifmag
2865 : : !! NAME
2866 : : !! calmaxdifmag
2867 : : !!
2868 : : !! FUNCTION
2869 : : !! Compute the maximum absolute value of the magnetization among all atoms,
2870 : : !! and the maximum absolute change in magnetization between the current and previous SCF cycles.
2871 : : !!
2872 : : !! INPUTS
2873 : : !! intgden(cplex,nspden,natom) = Integrated magnetic moments at the present SCF cycle.
2874 : : !! intgden0(cplex,nspden,natom) = Integrated magnetic moments at the previous SCF cycle.
2875 : : !! natom = Number of atoms.
2876 : : !! nspden = Number of spin-density components (typically 1, 2, or 4).
2877 : : !!
2878 : : !! OUTPUT
2879 : : !! maxmag = Maximum absolute value of magnetization among all atoms.
2880 : : !! difmag = Maximum absolute change in magnetization between present and previous SCF cycle.
2881 : : !!
2882 : : !! SIDE EFFECTS
2883 : : !! None
2884 : : !!
2885 : : !! NOTES
2886 : : !! This routine can handle both real and complex spin densities.
2887 : : !!
2888 : : !! SOURCE
2889 : :
2890 : 9241 : subroutine calmaxdifmag(cplex,intgden,intgden0,natom,nspden,maxmag,difmag)
2891 : :
2892 : : !Arguments ---------------------------------------------
2893 : : integer,intent(in) :: natom,nspden,cplex
2894 : : real(dp),intent(in) :: intgden(cplex,nspden,natom),intgden0(cplex,nspden,natom)
2895 : : real(dp),intent(out)::maxmag,difmag
2896 : : !Local variables ------------------------------
2897 : : integer :: iatom
2898 : : real(dp)::mag,mag0
2899 : 9241 : maxmag=zero;difmag=zero
2900 [ + + ]: 9241 : if (nspden==2 ) then
2901 [ + + ]: 19858 : do iatom=1,natom
2902 : 13312 : mag=intgden(1,1,iatom)-intgden(1,2,iatom)
2903 : 13312 : mag0=intgden0(1,1,iatom)-intgden0(1,2,iatom)
2904 : 13312 : maxmag=max(maxmag,abs(mag))
2905 : 13312 : difmag=max(difmag,abs(mag-mag0))
2906 [ + + ]: 19858 : if (cplex==2) then
2907 : 121 : mag=intgden(2,1,iatom)-intgden(2,2,iatom)
2908 : 121 : mag0=intgden(2,1,iatom)-intgden0(2,2,iatom)
2909 : 121 : maxmag=max(maxmag,abs(mag))
2910 : 121 : difmag=max(difmag,abs(mag-mag0))
2911 : : endif
2912 : : end do
2913 [ + - ]: 2695 : else if (nspden==4 ) then
2914 [ + + ]: 8520 : do iatom=1,natom
2915 [ + + ]: 23300 : mag =sqrt(sum(intgden(1,2:4,iatom)**2))
2916 [ + + ]: 23300 : mag0=sqrt(sum(intgden0(1,2:4,iatom)**2))
2917 : 5825 : maxmag=max(maxmag,abs(mag))
2918 : 5825 : difmag=max(difmag,abs(mag-mag0))
2919 [ + + ]: 8520 : if (cplex==2) then
2920 [ + + ]: 1240 : mag =sqrt(sum(intgden(2,2:4,iatom)**2))
2921 [ + + ]: 1240 : mag0=sqrt(sum(intgden0(2,2:4,iatom)**2))
2922 : 310 : maxmag=max(maxmag,abs(mag))
2923 : 310 : difmag=max(difmag,abs(mag-mag0))
2924 : : endif
2925 : : end do
2926 : : endif
2927 [ + + ]: 9241 : if (maxmag < tol8) then
2928 : 2027 : maxmag=0
2929 [ + + ]: 2027 : if (difmag < tol8) difmag=0
2930 : : endif
2931 : 9241 : end subroutine calmaxdifmag
2932 : : !!***
2933 : :
2934 : : !!****f* m_dens/fatsph_recip
2935 : : !! NAME
2936 : : !! fatsph_recip
2937 : : !!
2938 : : !! FUNCTION
2939 : : !! Compute the atomic spheres functions in reciprocal space as a product of
2940 : : !! a Bessel function times a Gaussian smearing for the boundary. The functions
2941 : : !! are subsequently Fourier transformed to real space.
2942 : : !!
2943 : : !! INPUTS
2944 : : !! mpi_enreg=information about MPI parallelization
2945 : : !! natom=number of atoms in cell.
2946 : : !! nfft=(effective) number of FFT grid points (for this processor)
2947 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2948 : : !! ntypat=number of atom types
2949 : : !! ratsm=smearing width for ratsph
2950 : : !! ratsph(ntypat)=radius of spheres around atoms
2951 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2952 : : !! typat(natom)=type of each atom
2953 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
2954 : : !!
2955 : : !! OUTPUT
2956 : : !! fatsph(nfft,natom)= functions defining the atomic spheres of integration in real space
2957 : : !! fatsphi3(n1,n2,n3,natom)= Same functions with triple indexing
2958 : : !!
2959 : : !! SOURCE
2960 : :
2961 : 0 : subroutine fatsph_recip(fatsph,fatsph3i,gmet,mpi_enreg,natom,nfft,ngfft,ntypat,&
2962 : 0 : & ratsm,ratsph,typat,ucvol,xred)
2963 : :
2964 : : !Arguments ------------------------------------
2965 : : !scalars
2966 : : integer,intent(in) :: natom,nfft,ntypat
2967 : : real(dp),intent(in) :: ratsm,ucvol
2968 : : type(MPI_type),intent(in) :: mpi_enreg
2969 : : !arrays
2970 : : integer,intent(in) :: ngfft(18),typat(natom)
2971 : : real(dp),intent(in) :: gmet(3,3),ratsph(ntypat)
2972 : : real(dp),intent(in) :: xred(3,natom)
2973 : : real(dp),intent(out):: fatsph(nfft,natom)
2974 : : real(dp),intent(out):: fatsph3i(ngfft(1),ngfft(2),ngfft(3),natom)
2975 : :
2976 : : !Local variables ------------------------------
2977 : : !scalars
2978 : : integer :: iatom
2979 : : integer :: i1,i2,i3,id1,id2,id3,ig1,ig2,ig3,ii,ii1,n1,n2,n3,n4,n5,n6
2980 : : real(dp) :: arg1,arg2,fac1,fac2,fac3,gq1,gq2,gq3,gcube
2981 : : real(dp) :: gsquar,gmag,gmagrad,rad,sfr,sfi,widthsq
2982 : : !arrays
2983 : 0 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
2984 : 0 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
2985 : : real(dp) :: gq(3)
2986 : 0 : real(dp) :: work1(2,nfft),work2(nfft,1),work3(ngfft(1),ngfft(2),ngfft(3),1)
2987 : :
2988 : : !******************************************************************
2989 : :
2990 : : !Geometric parameters
2991 : 0 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
2992 : 0 : n4=ngfft(4);n5=ngfft(5);n6=ngfft(6)
2993 : 0 : id1=n1/2+2
2994 : 0 : id2=n2/2+2
2995 : 0 : id3=n3/2+2
2996 : 0 : widthsq=ratsm**2
2997 : :
2998 : : !Get the distrib associated with this fft_grid
2999 : 0 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
3000 : :
3001 [ # # ]: 0 : do iatom=1, natom
3002 : :
3003 : 0 : ii=0
3004 [ # # # # ]: 0 : work1(:,:)=zero
3005 : : !G=0 term
3006 : 0 : rad=ratsph(typat(iatom))
3007 : 0 : work1(1,1)=four_pi*rad**3/(three*ucvol)
3008 [ # # ]: 0 : do i3=1,n3
3009 : 0 : ig3=i3-(i3/id3)*n3-1
3010 : 0 : gq3=dble(ig3)
3011 : 0 : gq(3)=gq3
3012 [ # # ]: 0 : do i2=1,n2
3013 [ # # ]: 0 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
3014 : 0 : ig2=i2-(i2/id2)*n2-1
3015 : 0 : gq2=dble(ig2)
3016 : 0 : gq(2)=gq2
3017 : :
3018 : : ! Note the lower limit of the next loop
3019 : 0 : ii1=1
3020 [ # # # # ]: 0 : if(i3==1 .and. i2==1 .and. ig2==0 .and. ig3==0)then
3021 : 0 : ii1=2
3022 : 0 : ii=ii+1
3023 : : end if
3024 [ # # ]: 0 : do i1=ii1,n1
3025 : 0 : ig1=i1-(i1/id1)*n1-1
3026 : 0 : gq1=dble(ig1)
3027 : 0 : gq(1)=gq1
3028 : 0 : ii=ii+1
3029 : :
3030 : 0 : gsquar=gsq_vl3(gq1,gq2,gq3)
3031 : 0 : gmag=sqrt(gsquar)
3032 : 0 : gcube=gmag*gsquar
3033 : 0 : gmagrad=two_pi*gmag*rad
3034 : 0 : arg1=-gsquar*pi**2*widthsq
3035 [ # # ]: 0 : arg2=two_pi*dot_product(xred(:,iatom),gq)
3036 : :
3037 : 0 : fac1=two/(two_pi**2*gcube*ucvol)
3038 : 0 : fac2=sin(gmagrad)-gmagrad*cos(gmagrad)
3039 : 0 : fac3=exp(arg1)
3040 : 0 : sfr=cos(arg2)
3041 : 0 : sfi=-sin(arg2)
3042 : :
3043 : 0 : work1(1,ii)=fac1*fac2*fac3*sfr
3044 : 0 : work1(2,ii)=fac1*fac2*fac3*sfi
3045 : :
3046 : : end do
3047 : : end if
3048 : : end do
3049 : : end do
3050 : :
3051 : : ! Transform to real space
3052 : 0 : call fourdp(1,work1,work2,1,mpi_enreg,nfft,1,ngfft,0)
3053 [ # # ]: 0 : fatsph(:,iatom)=work2(:,1)
3054 : :
3055 : 0 : call fftpac(1,mpi_enreg,1,n1,n2,n3,n1,n2,n3,ngfft,work2,work3,2)
3056 [ # # # # : 0 : fatsph3i(:,:,:,iatom)=work3(:,:,:,1)
# # ]
3057 : :
3058 : : end do !iatom
3059 : :
3060 : : contains
3061 : :
3062 : 0 : function gsq_vl3(g1,g2,g3)
3063 : :
3064 : : real(dp) :: gsq_vl3
3065 : : real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions
3066 : : !Define G^2 based on G space metric gmet.
3067 : : gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+&
3068 : : & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+&
3069 : 0 : & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1)
3070 : 0 : end function gsq_vl3
3071 : :
3072 : : end subroutine fatsph_recip
3073 : : !!***
3074 : :
3075 : : !!****f* m_dens/magmom_to_d2
3076 : : !! NAME
3077 : : !! magmom_to_d2
3078 : : !!
3079 : : !! FUNCTION
3080 : : !! Incorporates the magnetic moments in the ddb files as second
3081 : : !! order energy derivatives with respect to (ipert,idir) and a
3082 : : !! Zeeman field. Both total, i.e., response to a uniform Zeeman
3083 : : !! field (ipert=natom+5) and local (ipert=natom+11+1:2*natom+11)
3084 : : !! magnetic moments are considered. The Zeeman field directions
3085 : : !! are passed in Cartesian format.
3086 : : !! It also incorporates the second-order energy derivatives involving
3087 : : !! a scalar-potential perturbation from the charge-induced by
3088 : : !! another (ipert,idir). TODO: Check the consistency of signs in this case.
3089 : : !!
3090 : : !! INPUTS
3091 : : !! blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
3092 : : !! idir=direction of the perturbation
3093 : : !! intgden(cplex,nspden, natom)=integrated rhor or potential residual, for each atom in a sphere of radius ratsph.
3094 : : !! Representation differs according to nspden :
3095 : : !! if nspden=1, total density
3096 : : !! if nspden=2, spin up, then spin down
3097 : : !! if nspden=4, total density, then mag_x, mag_y, mag_z
3098 : : !! ipert=type of perturbation
3099 : : !! mpert=maximum number of perturbations
3100 : : !! natom=number of atoms in cell.
3101 : : !! nspden=number of spin-density components
3102 : : !! rhomag(2,nspden)=integral of charge or magnetization over the whole cell (also taking into account a possible imaginary part for DFPT).
3103 : : !!
3104 : : !! OUTPUT
3105 : : !! d2lo(2,3,mpert,3,mpert)= Local contributions to the second-order energy functional.
3106 : : !!
3107 : : !! SOURCE
3108 : :
3109 : 12 : subroutine magmom_to_d2(blkflg,cplex,d2lo,idir,intgden,ipert,mpert,natom,nspden,rhomag)
3110 : :
3111 : : !Arguments ---------------------------------------------
3112 : : !scalars
3113 : : integer,intent(in) :: cplex,idir,ipert,mpert,natom,nspden
3114 : : !arrays
3115 : : integer,intent(inout) :: blkflg(3,mpert,3,mpert)
3116 : : real(dp),intent(in) :: intgden(cplex,nspden,natom)
3117 : : real(dp),intent(in) :: rhomag(2,nspden)
3118 : : real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert)
3119 : : !Local variables ------------------------------
3120 : : !scalars
3121 : : integer :: iatom
3122 : :
3123 : : ! *************************************************************************
3124 : :
3125 : : ! We store in DDB the second-order energy derivatives, hence the negative
3126 : : ! sign applied to the induced magnetic moments.
3127 : :
3128 : : ! Incorporate total charge and magnetic moments
3129 [ - + ]: 12 : if (nspden==2) then
3130 : 0 : blkflg(1,natom+6,idir,ipert)= 1
3131 : 0 : d2lo(1,1,natom+6,idir,ipert)= -rhomag(1,1)
3132 [ # # ]: 0 : if (cplex==2) d2lo(2,1,natom+6,idir,ipert)= rhomag(2,1)
3133 : 0 : blkflg(3,natom+5,idir,ipert)= 1
3134 : 0 : d2lo(1,3,natom+5,idir,ipert)= -half*rhomag(1,2)
3135 [ # # ]: 0 : if (cplex==2) d2lo(2,3,natom+5,idir,ipert)= -half*rhomag(2,2)
3136 [ + - ]: 12 : else if (nspden==4) then
3137 : 12 : blkflg(1,natom+6,idir,ipert)=1
3138 : 12 : d2lo(1,1,natom+6,idir,ipert)= -rhomag(1,1)
3139 [ + + ]: 12 : if (cplex==2) d2lo(2,1,natom+6,idir,ipert)= rhomag(2,1)
3140 [ + + ]: 48 : blkflg(1:3,natom+5,idir,ipert)=1
3141 [ + + ]: 48 : d2lo(1,1:3,natom+5,idir,ipert)= -half*rhomag(1,2:4)
3142 [ + + + + ]: 36 : if (cplex==2) d2lo(2,1:3,natom+5,idir,ipert)= -half*rhomag(2,2:4)
3143 : : end if
3144 : :
3145 : : ! Incorporate local magnetic moments
3146 [ - + ]: 12 : if (nspden==2) then
3147 [ # # ]: 0 : do iatom= 1, natom
3148 : 0 : blkflg(3,natom+11+iatom,idir,ipert)= 1
3149 : 0 : d2lo(1,3,natom+11+iatom,idir,ipert)= -half*intgden(1,2,iatom)
3150 [ # # ]: 0 : if (cplex==2) d2lo(2,3,natom+11+iatom,idir,ipert)= -half*intgden(2,2,iatom)
3151 : : end do
3152 [ + - ]: 12 : else if (nspden==4) then
3153 [ + + ]: 84 : do iatom= 1, natom
3154 [ + + ]: 288 : blkflg(1:3,natom+11+iatom,idir,ipert)= 1
3155 [ + + ]: 288 : d2lo(1,1:3,natom+11+iatom,idir,ipert)= -half*intgden(1,2:4,iatom)
3156 [ + + + + ]: 228 : if (cplex==2) d2lo(2,1:3,natom+11+iatom,idir,ipert)= -half*intgden(2,2:4,iatom)
3157 : : end do
3158 : : end if
3159 : :
3160 : 12 : end subroutine magmom_to_d2
3161 : : !!***
3162 : :
3163 [ # # # # : 0 : end module m_dens
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
3164 : : !!***
|
