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 1413 : subroutine add_atomic_fcts(natom,nspden,rprimd,mpi_enreg,nfft,ngfft,ntypat,option,ratsph, &
467 1413 : 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 1413 : real(dp), allocatable :: difx(:)
497 : real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3)
498 : real(dp) :: tsec(2)
499 1413 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
500 1413 : 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 1413 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
505 :
506 1413 : n1 = ngfft(1)
507 1413 : n2 = ngfft(2)
508 1413 : 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 1413 : if(option==0)then
514 1425865 : nv_constr_dft_r = zero
515 : endif
516 :
517 : !Get the distrib associated with this fft_grid
518 1413 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
519 :
520 : !Loop over atoms
521 : !-------------------------------------------
522 4664 : do iatom=1,natom
523 :
524 14517 : if(sum(coeffs_constr_dft(1:nspden,iatom)**2)<tol12)then
525 : cycle
526 : endif
527 :
528 : ! Define a "box" around the atom
529 3025 : r2atsph=1.0000001_dp*ratsph(typat(iatom))**2
530 3025 : rr1=sqrt(r2atsph*gmet(1,1))
531 3025 : rr2=sqrt(r2atsph*gmet(2,2))
532 3025 : rr3=sqrt(r2atsph*gmet(3,3))
533 :
534 3025 : n1a=int((xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1
535 3025 : n1b=int((xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1
536 3025 : n2a=int((xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2
537 3025 : n2b=int((xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2
538 3025 : n3a=int((xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3
539 3025 : n3b=int((xred(3,iatom)+rr3+ishift)*n3 )-ishift*n3
540 :
541 3025 : ratsm2 = -(ratsm**2 - 2*ratsph(typat(iatom))*ratsm)
542 :
543 9075 : ABI_MALLOC(difx,(n1a:n1b))
544 30554 : do i1=n1a,n1b
545 30554 : difx(i1)=dble(i1)/dble(n1)-xred(1,iatom)
546 : enddo ! i1
547 :
548 29416 : do i3=n3a,n3b
549 26391 : iz=mod(i3+ishift*n3,n3)
550 29416 : if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
551 26391 : izloc = ffti3_local(iz+1) - 1
552 26391 : difz=dble(i3)/dble(n3)-xred(3,iatom)
553 383456 : do i2=n2a,n2b
554 357065 : iy=mod(i2+ishift*n2,n2)
555 357065 : dify=dble(i2)/dble(n2)-xred(2,iatom)
556 357065 : rx23=dify*rprimd(1,2)+difz*rprimd(1,3)
557 357065 : ry23=dify*rprimd(2,2)+difz*rprimd(2,3)
558 357065 : rz23=dify*rprimd(3,2)+difz*rprimd(3,3)
559 357065 : r2_23=rx23**2+ry23**2+rz23**2
560 357065 : r2_11=rprimd(1,1)**2+rprimd(2,1)**2+rprimd(3,1)**2
561 357065 : r2_123=2*(rprimd(1,1)*rx23+rprimd(2,1)*ry23+rprimd(3,1)*rz23)
562 6713401 : do i1=n1a,n1b
563 6329945 : r2=(difx(i1)*r2_11+r2_123)*difx(i1)+r2_23
564 6329945 : if (r2 > r2atsph) cycle
565 3369523 : call radsmear(dfsm,fsm,r2,r2atsph,ratsm2)
566 3369523 : ix=mod(i1+ishift*n1,n1)
567 : ! Identify the fft indexes of the rectangular grid around the atom
568 3369523 : ifft_local=1+ix+n1*(iy+n2*izloc)
569 11585447 : 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 4664 : 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 1413 : 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 1413 : 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 6923 : subroutine constrained_dft_free(constrained_dft)
738 :
739 : !Arguments ------------------------------------
740 : type(constrained_dft_t),intent(inout):: constrained_dft
741 : ! ***********************************************************************************************
742 :
743 6923 : ABI_SFREE(constrained_dft%chrgat)
744 6923 : ABI_SFREE(constrained_dft%constraint_kind)
745 6923 : ABI_SFREE(constrained_dft%intgf2)
746 6923 : ABI_SFREE(constrained_dft%ratsph)
747 6923 : ABI_SFREE(constrained_dft%spinat)
748 6923 : ABI_SFREE(constrained_dft%typat)
749 6923 : ABI_SFREE(constrained_dft%ziontypat)
750 6923 : ABI_SFREE(constrained_dft%znucl)
751 :
752 6923 : 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 1304 : 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 2608 : integer :: ipiv(c_dft%natom)
828 : real(dp) :: corr_denmag(4),gr_intgd(3),strs_intgd(6)
829 1304 : real(dp), allocatable :: coeffs_constr_dft(:,:) ! nspden,natom
830 1304 : real(dp), allocatable :: gr_intgden(:,:,:) ! 3,nspden,natom
831 1304 : real(dp), allocatable :: intgden(:,:) ! nspden,natom
832 1304 : real(dp), allocatable :: intgden_delta(:,:) ! nspden,natom
833 1304 : real(dp), allocatable :: intgres_tmp(:,:) ! nspden,natom
834 1304 : real(dp), allocatable :: intgr(:,:) ! natom,nspden
835 1304 : real(dp), allocatable :: strs_intgden(:,:,:) ! 6,nspden,natom
836 2608 : 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 1304 : natom=c_dft%natom
846 1304 : nfftf=c_dft%nfftf
847 1304 : nspden=c_dft%nspden
848 1304 : ntypat=c_dft%ntypat
849 :
850 : !We need the integrated magnetic moments
851 5216 : ABI_MALLOC(intgden,(nspden,natom))
852 5216 : ABI_MALLOC(gr_intgden,(3,nspden,natom))
853 5216 : 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 1304 : xred,1,cplex1,qgbt,use_gbt,intgden=intgden,gr_intgden=gr_intgden,rhomag=rhomag,strs_intgden=strs_intgden)
857 2608 : 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 3912 : ABI_MALLOC(intgres_tmp,(nspden,natom))
865 15364 : intgres_tmp(:,:)=zero
866 : call calcdenmagsph(mpi_enreg,natom,nfftf,c_dft%ngfftf,nspden,ntypat,&
867 1304 : & 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 15364 : intgres(:,:)=zero
875 4446 : do iatom=1,natom
876 3142 : if(nspden==1)then
877 210 : intgres(1,iatom)=intgres_tmp(1,iatom)
878 2932 : 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 2422 : else if(nspden==4)then
882 : !Change the potential residual to the density+magnetization convention
883 2422 : intgres(1,iatom)=half*(intgres_tmp(1,iatom)+intgres_tmp(2,iatom))
884 2422 : intgres(2,iatom)= intgres_tmp(3,iatom)
885 2422 : intgres(3,iatom)=-intgres_tmp(4,iatom)
886 2422 : intgres(4,iatom)=half*(intgres_tmp(1,iatom)-intgres_tmp(2,iatom))
887 : endif
888 3142 : conkind=c_dft%constraint_kind(c_dft%typat(iatom))
889 3142 : if(conkind <10)intgres(1,iatom)=zero
890 4671 : if( mod(conkind,10)==0 .and. nspden>1)intgres(2:nspden,iatom)=zero
891 : enddo
892 : !Print the potential residuals
893 2608 : 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 1304 : ABI_FREE(intgres_tmp)
895 :
896 : !Also exchanges the spin and atom indices to prepare the solution of the linear system of equation
897 3912 : ABI_MALLOC(intgr,(natom,nspden))
898 4446 : do iatom=1,natom
899 15364 : 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 3912 : do ii=1,2 ! Charge, then spin
907 25448 : intgf2(:,:)=zero
908 8892 : do iatom=1,natom
909 8892 : intgf2(iatom,iatom)=c_dft%intgf2(iatom,iatom)
910 : enddo
911 6284 : do iatom=1,natom-1
912 15096 : do jatom=iatom,natom
913 12488 : if(c_dft%intgf2(iatom,jatom)>tol8)then
914 : !In the charge case, must have both atoms with constraints bigger than 10
915 5900 : if(ii==1)then
916 2950 : 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 5900 : if(ii==2)then
924 2950 : 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 2738 : intgf2(iatom,jatom)=c_dft%intgf2(iatom,jatom)
927 2738 : 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 2608 : call dsytrf('U',natom,intgf2,natom,ipiv,work,2*natom,info)
936 : ! call dsytri('U',natom,intgf2,natom,ipiv,work,info)
937 2608 : if(ii==1)then
938 1304 : call dsytrs('U',natom,1,intgf2,natom,ipiv,intgr,natom,info)
939 1304 : else if(ii==2 .and. nspden>1)then
940 1199 : 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 10196 : do iatom=1,natom
945 6284 : if(ii==1)intgres(1,iatom)=intgr(iatom,1)
946 16668 : 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 1304 : e_constrained_dft=zero
960 13872 : grcondft=zero
961 1304 : strscondft=zero
962 5216 : ABI_MALLOC(intgden_delta,(nspden,natom))
963 15364 : intgden_delta(:,:)=zero
964 4446 : do iatom=1,natom
965 :
966 : ! The integrated density must be in the total density+magnetization representation
967 3142 : 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 3142 : conkind=c_dft%constraint_kind(c_dft%typat(iatom))
985 :
986 3142 : 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 3142 : if( mod(conkind,10)==1 .and. nspden>1)then
995 : !Fix the different components of the magnetization vector
996 2458 : if(nspden==2)intgden_delta(2,iatom)=intgden(2,iatom)-c_dft%spinat(3,iatom)
997 8908 : 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 14060 : e_constrained_dft=e_constrained_dft-sum(intgden_delta(:,iatom)*intgres(:,iatom))
1073 12568 : do ii=1,3
1074 45322 : 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 23298 : do ii=1,6
1078 87502 : 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 1304 : ABI_FREE(gr_intgden)
1093 1304 : ABI_FREE(strs_intgden)
1094 :
1095 3912 : ABI_MALLOC(coeffs_constr_dft,(nspden,natom))
1096 15364 : 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 4446 : do iatom=1,natom
1101 :
1102 : !Computation of the correction in terms of density and magnetization coefficients.
1103 3142 : conkind=c_dft%constraint_kind(c_dft%typat(iatom))
1104 3142 : corr_denmag(:)=zero
1105 :
1106 3142 : if(conkind >=10)then
1107 825 : corr_denmag(1)=intgden_delta(1,iatom)*c_dft%magcon_lambda - intgres(1,iatom)
1108 : endif
1109 :
1110 3142 : if( mod(conkind,10)==1 .and. nspden>1)then
1111 :
1112 : !Fix the different components of the magnetization vector
1113 2458 : if(nspden==2)corr_denmag(2)=intgden_delta(2,iatom)*c_dft%magcon_lambda - intgres(2,iatom)
1114 8908 : 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 4446 : if(nspden==1)then
1150 : !From charge to potential
1151 210 : coeffs_constr_dft(1,iatom)=corr_denmag(1)
1152 2932 : else if(nspden==2 .or. nspden==4)then
1153 : !From charge and magnetization to potential
1154 2932 : coeffs_constr_dft(1,iatom)=corr_denmag(1)+corr_denmag(nspden)
1155 2932 : coeffs_constr_dft(2,iatom)=corr_denmag(1)-corr_denmag(nspden)
1156 2932 : if(nspden==4)then
1157 2422 : coeffs_constr_dft(3,iatom)= corr_denmag(2)
1158 2422 : 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 1304 : option=1
1166 : call add_atomic_fcts(natom,nspden,c_dft%rprimd,mpi_enreg,nfftf,c_dft%ngfftf,ntypat,option,&
1167 1304 : & c_dft%ratsph,c_dft%ratsm,c_dft%typat,coeffs_constr_dft,vresid,xred) ! Ben change: add ratsm
1168 :
1169 1304 : ABI_FREE(coeffs_constr_dft)
1170 1304 : ABI_FREE(intgden)
1171 1304 : ABI_FREE(intgden_delta)
1172 1304 : ABI_FREE(intgr)
1173 :
1174 1304 : 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 56610 : subroutine calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsm,ratsph,rhor,rprimd,typat,xred,&
1530 21602 : & 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 37740 : 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 18870 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1569 37740 : integer :: overlap_ij(natom,natom)
1570 : real(dp) :: gmet(3,3),gprimd(3,3),gr_intg(3,4)
1571 37740 : 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 56610 : real(dp) :: dist_ij(natom,natom),intgden_(cplex,nspden,natom)!,intgden_im_(nspden,natom)
1575 37740 : real(dp) :: my_xred(3, natom), rmet(3,3),xshift(3, natom), taumr_local(3)
1576 18870 : real(dp), allocatable :: fsm_atom(:,:)
1577 18870 : real(dp), ABI_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 18870 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
1584 18870 : n4=ngfft(4);n5=ngfft(5);n6=ngfft(6)
1585 18870 : nfftot=n1*n2*n3
1586 :
1587 : !Manage optinal arguments
1588 258634 : if(present(intgden)) intgden=zero
1589 65684 : if(present(gr_intgden)) gr_intgden=zero
1590 98438 : if(present(strs_intgden)) strs_intgden=zero
1591 18870 : qphon_=zero
1592 18870 : if(present(qphon))then
1593 74 : qphon_=qphon
1594 : endif
1595 :
1596 18870 : if(present(fatsph)) then
1597 62 : fatsph_ => fatsph
1598 952754 : fatsph_=zero
1599 : end if
1600 18870 : if(present(taumr)) then
1601 62 : taumr_ => taumr
1602 2858324 : taumr_=zero
1603 : end if
1604 :
1605 18870 : 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 203078 : 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 18870 : 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 18870 : grid_found=.false.
1637 :
1638 18870 : if(n2 == mpi_enreg%distribfft%n2_coarse ) then
1639 13597 : if(n3== size(mpi_enreg%distribfft%tab_fftdp3_distrib)) then
1640 13543 : fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3_distrib
1641 13543 : ffti3_local => mpi_enreg%distribfft%tab_fftdp3_local
1642 13543 : grid_found=.true.
1643 : end if
1644 : end if
1645 :
1646 18870 : if(n2 == mpi_enreg%distribfft%n2_fine ) then
1647 5535 : if(n3 == size(mpi_enreg%distribfft%tab_fftdp3dg_distrib)) then
1648 5535 : fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3dg_distrib
1649 5535 : ffti3_local => mpi_enreg%distribfft%tab_fftdp3dg_local
1650 : grid_found = .true.
1651 : end if
1652 : end if
1653 :
1654 13335 : ABI_CHECK(grid_found, "Unable to find an allocated distrib for this fft grid")
1655 :
1656 : !Loop over atoms
1657 : !-------------------------------------------
1658 64922 : do iatom=1,natom
1659 :
1660 : ! Define a "box" around the atom that extends outside of the unit cell
1661 46052 : r2atsph=1.0000001_dp*ratsph(typat(iatom))**2
1662 46052 : rr1=sqrt(r2atsph*gmet(1,1))
1663 46052 : rr2=sqrt(r2atsph*gmet(2,2))
1664 46052 : rr3=sqrt(r2atsph*gmet(3,3))
1665 :
1666 46052 : n1a=int((my_xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1
1667 46052 : n1b=int((my_xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1
1668 46052 : n2a=int((my_xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2
1669 46052 : n2b=int((my_xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2
1670 46052 : n3a=int((my_xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3
1671 46052 : 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 46052 : ratsm2 = (2*ratsph(typat(iatom))-ratsm)*ratsm
1675 :
1676 417416 : intg(:,:)=zero
1677 46052 : gr_intg(:,:)=zero
1678 46052 : strs_intg(:,:)=zero
1679 :
1680 629449 : do i3=n3a,n3b
1681 583397 : iz=mod(i3+ishift*n3,n3)
1682 :
1683 629449 : if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
1684 :
1685 572411 : izloc = ffti3_local(iz+1) - 1
1686 572411 : difz=dble(i3)/dble(n3)-my_xred(3,iatom)
1687 9884125 : do i2=n2a,n2b
1688 9311714 : iy=mod(i2+ishift*n2,n2)
1689 9311714 : dify=dble(i2)/dble(n2)-my_xred(2,iatom)
1690 208950021 : do i1=n1a,n1b
1691 199065896 : ix=mod(i1+ishift*n1,n1)
1692 199065896 : 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 199065896 : 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 199065896 : ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1708 199065896 : rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1709 199065896 : r2=rx**2+ry**2+rz**2
1710 :
1711 : ! Identify the fft indexes of the rectangular grid around the atom
1712 199065896 : if(r2 > r2atsph) then
1713 : cycle
1714 : end if
1715 :
1716 88854453 : call radsmear(dfsm,fsm,r2,r2atsph,ratsm2)
1717 :
1718 88854453 : ifft_local=1+ix+n1*(iy+n2*izloc)
1719 :
1720 88854453 : if(present(fatsph)) then
1721 60636 : fatsph_(ifft_local,iatom)=fsm
1722 : end if
1723 88854453 : 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 88854453 : taumr_local(1)=difx
1731 88854453 : taumr_local(2)=dify
1732 88854453 : taumr_local(3)=difz
1733 355417812 : arg=two_pi*dot_product(qphon_,taumr_local)
1734 88854453 : phr1d_re=dcos(arg)
1735 88854453 : phr1d_im=dsin(arg)
1736 88854453 : ifft_local_cplex=1+cplex*(ifft_local-1)
1737 :
1738 88854453 : 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 88854453 : if (nspden==4 .and. use_gbt /= 0) then
1745 : ! If GBT is on, we have to include the e^{iq.r} phase
1746 206397 : qr = two_pi * (qgbt(1)*dble(i1)/dble(n1) + qgbt(2)*dble(i2)/dble(n2) + qgbt(3)*dble(i3)/dble(n3))
1747 206397 : rhor_local(1) = rhor(ifft_local,1)
1748 206397 : mx = rhor(ifft_local,2)
1749 206397 : my = rhor(ifft_local,3)
1750 206397 : mz = rhor(ifft_local,4)
1751 : ! e^-iqr :cos(qr)-isin(qr)
1752 206397 : rhor_local(2) = cos(qr)*mx - sin(qr)*my ! mx
1753 206397 : rhor_local(3) = sin(qr)*mx + cos(qr)*my ! my
1754 206397 : rhor_local(4) = mz
1755 : else
1756 259685971 : rhor_local(1:nspden) = rhor(ifft_local,1:nspden)
1757 : end if
1758 : ! Integral of density or potential residual
1759 88854453 : if (cplex==1) then
1760 258756197 : intg(1,1:nspden)=intg(1,1:nspden)+fsm*rhor(ifft_local,1:nspden)
1761 457533 : else if (cplex==2) then
1762 1830132 : if (sum(qphon_(:)**2)<tol8) then
1763 1823959 : intg(1,1:nspden)=intg(1,1:nspden)+fsm*rhor(2*ifft_local-1,1:nspden)
1764 1823959 : intg(2,1:nspden)=intg(2,1:nspden)+fsm*rhor(2*ifft_local ,1:nspden)
1765 : else
1766 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)
1767 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)
1768 : end if
1769 : end if
1770 187020620 : if((present(gr_intgden).or.present(strs_intgden)).and. option<10 .and. ratsm2>tol12)then
1771 11337976 : do ispden=1,nspden
1772 7873026 : fact=dfsm*rhor(ifft_local_cplex,ispden)
1773 7873026 : if(present(gr_intgden))then
1774 7873026 : gr_intg(1,ispden)=gr_intg(1,ispden)+difx*fact
1775 7873026 : gr_intg(2,ispden)=gr_intg(2,ispden)+dify*fact
1776 7873026 : gr_intg(3,ispden)=gr_intg(3,ispden)+difz*fact
1777 : endif
1778 11337976 : if(present(strs_intgden))then
1779 7873026 : strs_intg(1,ispden)=strs_intg(1,ispden)+difx*difx*fact
1780 7873026 : strs_intg(2,ispden)=strs_intg(2,ispden)+dify*dify*fact
1781 7873026 : strs_intg(3,ispden)=strs_intg(3,ispden)+difz*difz*fact
1782 7873026 : strs_intg(4,ispden)=strs_intg(4,ispden)+dify*difz*fact
1783 7873026 : strs_intg(5,ispden)=strs_intg(5,ispden)+difx*difz*fact
1784 7873026 : strs_intg(6,ispden)=strs_intg(6,ispden)+difx*dify*fact
1785 : endif
1786 : enddo
1787 : endif
1788 : end do
1789 : end do
1790 : end if
1791 : end do
1792 :
1793 :
1794 : !DEBUG
1795 : ! n1c=(n1b+n1a)/2
1796 : ! n2c=(n2b+n2a)/2
1797 : ! do i3= n3a-5,n3b+5
1798 : ! n1c=mod(n1c+ishift*n1,n1)
1799 : ! n2c=mod(n2c+ishift*n2,n2)
1800 : ! iz=mod(i3+ishift*n3,n3)
1801 : ! end do
1802 : !ENDDEBUG
1803 :
1804 :
1805 46052 : if(present(intgf2) .and. neighbor_overlap==0)then
1806 76 : intgf2(iatom,iatom)=intgf2(iatom,iatom)*ucvol/dble(nfftot)
1807 : endif
1808 :
1809 417416 : intg(:,:)=intg(:,:)*ucvol/dble(nfftot)
1810 :
1811 46052 : if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
1812 : ! Convert to gradient in reduced coordinates
1813 166526 : gr_intg=matmul(rmet,gr_intg)
1814 53414 : gr_intg(:,:)=-gr_intg(:,:)*two*ucvol/dble(nfftot)
1815 : !DEBUG
1816 : ! write(6,*)' calcdenmagsph : intg(1)=',intg(1)
1817 : ! write(6,*)' calcdenmagsph : iatom,gr_intg(1,1)=',iatom,gr_intg(1,1)
1818 : ! call flush(6)
1819 : ! stop
1820 : !ENDDEBUG
1821 : endif
1822 :
1823 46052 : if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then
1824 : ! Convert to stress in cartesian coordinates, for each spin constraint
1825 15710 : do isp=1,4
1826 : ! First change the representation
1827 12568 : strs(1,1)=strs_intg(1,isp) ; strs(2,2)=strs_intg(2,isp) ; strs(3,3)=strs_intg(3,isp)
1828 12568 : strs(2,3)=strs_intg(4,isp) ; strs(3,2)=strs_intg(4,isp)
1829 12568 : strs(1,3)=strs_intg(5,isp) ; strs(3,1)=strs_intg(5,isp)
1830 12568 : strs(1,2)=strs_intg(6,isp) ; strs(2,1)=strs_intg(6,isp)
1831 :
1832 : ! Then perform representation change, following Eq.(25) in Hamann2005
1833 50272 : do ii=1,3
1834 615832 : strs_cartred(:,ii)=matmul(rprimd,strs(:,ii))
1835 : enddo
1836 50272 : do ii=1,3
1837 615832 : strs(ii,:)=matmul(rprimd,strs_cartred(ii,:))
1838 : enddo
1839 :
1840 : !DEBUG
1841 : ! rprimd_mod=rprimd
1842 : ! rprimd_mod(1,1)=rprimd(1,1)*(one+strain)
1843 : ! do ii=1,3
1844 : ! strs_cartred(:,ii)=matmul(rprimd_mod,strs(:,ii))
1845 : ! enddo
1846 : ! do ii=1,3
1847 : ! strs(ii,:)=matmul(rprimd_mod,strs_cartred(ii,:))
1848 : ! enddo
1849 : !ENDDEBUG
1850 :
1851 12568 : strs_intg(1,isp)=two*strs(1,1) ; strs_intg(2,isp)=two*strs(2,2) ; strs_intg(3,isp)=two*strs(3,3)
1852 12568 : strs_intg(4,isp)=strs(2,3)+strs(3,2)
1853 12568 : strs_intg(5,isp)=strs(1,3)+strs(3,1)
1854 15710 : strs_intg(6,isp)=strs(1,2)+strs(2,1)
1855 : enddo
1856 91118 : strs_intg(:,:)=strs_intg(:,:)/dble(nfftot)
1857 : !DEBUG
1858 : ! strs_intg(:,:)=zero
1859 : ! write(6,*)' calcdenmagsph : iatom,-strs_intg(1,1)*ucvol=',iatom,strs_intg(1,1)*ucvol
1860 : ! call flush(6)
1861 : ! stop
1862 : !ENDDEBUG
1863 : endif
1864 :
1865 64922 : if(nspden==2 .and. option/=11)then
1866 : ! Specific treatment of collinear density, due to the storage mode.
1867 : ! intgden_(1,iatom)= integral of up density
1868 : ! intgden_(2,iatom)= integral of dn density
1869 29771 : intgden_(1:cplex,1,iatom)=intg(1:cplex,2)
1870 29771 : intgden_(1:cplex,2,iatom)=intg(1:cplex,1)-intg(1:cplex,2)
1871 : !if (cplex==2) then
1872 : ! intgden_im_(1,iatom)=intg(2,2)
1873 : ! intgden_im_(2,iatom)=intg(2,1)-intg(2,2)
1874 : !endif
1875 14825 : if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
1876 2040 : gr_intgden(:,1,iatom)=gr_intg(:,2)
1877 2040 : gr_intgden(:,2,iatom)=gr_intg(:,1)-gr_intg(:,2)
1878 : endif
1879 14825 : if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then
1880 3570 : strs_intgden(:,1,iatom)=strs_intg(:,2)
1881 3570 : strs_intgden(:,2,iatom)=strs_intg(:,1)-strs_intg(:,2)
1882 : endif
1883 : else
1884 166429 : intgden_(1:cplex,1:nspden,iatom)=intg(1:cplex,1:nspden)
1885 31227 : if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then
1886 42224 : gr_intgden(:,1:nspden,iatom)=gr_intg(:,1:nspden)
1887 : endif
1888 31227 : if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then
1889 71918 : strs_intgden(:,1:nspden,iatom)=strs_intg(:,1:nspden)
1890 : endif
1891 : endif
1892 :
1893 : end do ! iatom
1894 :
1895 : !-------------------------------------------
1896 : !
1897 : ! In case intgf2 must be computed, while the atoms overlap, a double loop over atoms is needed
1898 18870 : if(present(intgf2) .and. neighbor_overlap==1)then
1899 69 : do iatom=1,natom
1900 149 : do jatom=iatom,natom
1901 149 : if(overlap_ij(iatom,jatom)/=0)then
1902 1935 : do i3=1,n3
1903 1836 : iz=mod(i3,n3)
1904 1935 : if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
1905 1836 : izloc = ffti3_local(iz+1) - 1
1906 42300 : do i2=1,n2
1907 40464 : iy=mod(i2,n2)
1908 1128348 : do i1=1,n1
1909 1086048 : ix=mod(i1,n1)
1910 1086048 : ifft_local=1+ix+n1*(iy+n2*izloc)
1911 1126512 : intgf2(iatom,jatom)=intgf2(iatom,jatom)+fsm_atom(ifft_local,iatom)*fsm_atom(ifft_local,jatom)
1912 : enddo
1913 : enddo
1914 : endif
1915 : enddo ! i3
1916 99 : intgf2(iatom,jatom)=intgf2(iatom,jatom)*ucvol/dble(nfftot)
1917 : endif
1918 : enddo
1919 69 : if(iatom/=1)then
1920 80 : do jatom=1,iatom-1
1921 80 : intgf2(iatom,jatom)=intgf2(jatom,iatom)
1922 : enddo
1923 : endif
1924 : enddo
1925 19 : ABI_FREE(fsm_atom)
1926 : endif
1927 :
1928 : !MPI parallelization
1929 18870 : if(present(intgf2)) then
1930 54 : if(mpi_enreg%nproc_fft>1)then
1931 0 : call timab(48,1,tsec)
1932 0 : call xmpi_sum(intgf2,mpi_enreg%comm_fft,ierr)
1933 0 : call timab(48,2,tsec)
1934 : end if
1935 : end if
1936 :
1937 : !-------------------------------------------
1938 :
1939 : !MPI parallelization
1940 18870 : if(present(intgden) .or. option/=0) then
1941 18816 : if(mpi_enreg%nproc_fft>1)then
1942 544 : call timab(48,1,tsec)
1943 544 : call xmpi_sum(intgden_,mpi_enreg%comm_fft,ierr)
1944 544 : call timab(48,2,tsec)
1945 : end if
1946 258634 : if(present(intgden))intgden = intgden_
1947 : end if
1948 :
1949 18870 : if(present(gr_intgden) .and. option<10 .and. ratsm2>tol12) then
1950 1304 : if(mpi_enreg%nproc_fft>1)then
1951 0 : call timab(48,1,tsec)
1952 0 : call xmpi_sum(gr_intgden,mpi_enreg%comm_fft,ierr)
1953 0 : call timab(48,2,tsec)
1954 : end if
1955 : end if
1956 :
1957 18870 : if(present(strs_intgden) .and. option<10 .and. ratsm2>tol12) then
1958 1304 : if(mpi_enreg%nproc_fft>1)then
1959 0 : call timab(48,1,tsec)
1960 0 : call xmpi_sum(strs_intgden,mpi_enreg%comm_fft,ierr)
1961 0 : call timab(48,2,tsec)
1962 : end if
1963 : end if
1964 :
1965 : !EB - Compute magnetization of the whole cell
1966 18870 : if(present(dentot) .or. present(rhomag))then
1967 74090 : rhomag_(:,:)=zero
1968 11282 : if(nspden==2) then
1969 50425001 : do ifft=1,nfft
1970 50422382 : jfft=1+cplex*(ifft-1)
1971 101245604 : rhomag_(1:cplex,1)=rhomag_(1:cplex,1)+rhor(jfft:jfft+cplex-1,1) ! real & imag part of density
1972 101248223 : 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
1973 : end do
1974 8663 : else if(nspden==4) then
1975 27223670 : do ifft=1,nfft
1976 27221325 : jfft=1+cplex*(ifft-1)
1977 247890142 : rhomag_(1:cplex,1:nspden)=rhomag_(1:cplex,1:nspden)+rhor(jfft:jfft+cplex-1,1:nspden)
1978 : end do
1979 : end if
1980 :
1981 53872 : rhomag_(1:cplex,1:nspden)=rhomag_(1:cplex,1:nspden)*ucvol/dble(nfftot)
1982 :
1983 : !MPI parallelization
1984 11282 : if(mpi_enreg%nproc_fft>1)then
1985 440 : call timab(48,1,tsec)
1986 440 : call xmpi_sum(rhomag_,mpi_enreg%comm_fft,ierr)
1987 440 : call timab(48,2,tsec)
1988 : end if
1989 :
1990 14812 : if(present(dentot)) dentot(:)=rhomag_(1,:)
1991 74090 : if(present(rhomag)) rhomag(:,:)=rhomag_(:,:)
1992 : endif
1993 :
1994 : !DEBUG BUT KEEP
1995 : !if(.false.) call printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,'DEN.vtk')
1996 :
1997 59342 : end subroutine calcdenmagsph
1998 : !!***
1999 :
2000 : !!****f* m_dens/prtdenmagsph
2001 : !! NAME
2002 : !! prtdenmagsph
2003 : !!
2004 : !! FUNCTION
2005 : !! Print integral of total density inside spheres around atoms,
2006 : !! and optionally integral of potential residual (also gradient of the constraint energy wrt constraint).
2007 : !!
2008 : !! INPUTS
2009 : !! intgden(nspden, natom)=integrated rhor or potential residual, for each atom in a sphere of radius ratsph.
2010 : !! if option <10, intgden is a density (+magnetization)
2011 : !! Representation differs according to nspden :
2012 : !! if nspden=1, total density
2013 : !! if nspden=2, spin up, then spin down
2014 : !! if nspden=4, total density, then mag_x, mag_y, mag_z
2015 : !! if 20>option>=10, intgden is a potential residual (+spin magnetic field residual).
2016 : !! Representation is: first, mean potential; then (if nspden>=2) B_z for nspden=2, B_x, B_y and B_z for nspden=4.
2017 : !! if option>=20, intgden is a a gradient wrt target (=torque), also potential residual (+spin magnetic field residual) multiplied by f-1.
2018 : !! natom=number of atoms in cell.
2019 : !! nspden=number of spin-density components
2020 : !! ntypat=number of atom types
2021 : !! units=unit numbers for printing
2022 : !! 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..
2023 : !! When 1, 11, 21, the default printing is on (to unit nunit), if -1, 2, 3, 4, special printing options, if 0 no printing.
2024 : !! ratsm=smearing width for ratsph
2025 : !! ratsph(ntypat)=radius of spheres around atoms
2026 : !! rhomag(2,nspden)=integral of charge or magnetization over the whole cell (also taking into account a possible imaginary part for DFPT).
2027 : !! typat(natom)=type of each atom
2028 : !! ziontypat(ntypat)= --optional-- ionic charge of each atomic type
2029 : !!
2030 : !! OUTPUT
2031 : !! Printing
2032 : !!
2033 : !! SOURCE
2034 :
2035 9802 : subroutine prtdenmagsph(cplex, intgden, natom, nspden, ntypat, units, option, qgbt, ratsm, ratsph, rhomag, typat, znucl, spinaxis, ziontypat)
2036 :
2037 : !Arguments ---------------------------------------------
2038 : !scalars
2039 : integer,intent(in) :: natom,nspden,ntypat
2040 : integer,intent(in) :: units(:)
2041 : real(dp),intent(in) :: ratsm
2042 : integer ,intent(in) :: option
2043 : integer, intent(in) :: cplex
2044 : !arrays
2045 : integer,intent(in) :: typat(natom)
2046 : real(dp),intent(in) :: intgden(cplex,nspden,natom),qgbt(3),znucl(ntypat),spinaxis(3)
2047 : real(dp),intent(in) :: ratsph(ntypat),rhomag(2,nspden)
2048 : real(dp),intent(in),optional :: ziontypat(ntypat)
2049 :
2050 : !Local variables ------------------------------
2051 : !scalars
2052 : integer :: iatom,icplex
2053 : real(dp) :: mag_coll , mag_x, mag_y, mag_z ! EB
2054 : real(dp) :: mag_coll_im, mag_x_im, mag_y_im, mag_z_im ! SPr
2055 : real(dp) :: rho_tot, rho_tot_im
2056 : real(dp) :: sum_mag, sum_mag_x,sum_mag_y,sum_mag_z,sum_rho_up,sum_rho_dn,sum_rho_tot ! EB
2057 : 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
2058 : real(dp) :: mag_r, mag_theta, mag_phi, vec(3)
2059 : real(dp) :: sum_mag_r, sum_mag_theta, sum_mag_phi
2060 : real(dp) :: exact_mag_r, exact_mag_theta, exact_mag_phi
2061 : real(dp) :: exact_mag_r_im, exact_mag_theta_im, exact_mag_phi_im
2062 : 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)
2063 : character(len=500) :: msg,msg1
2064 : character(len=500) :: msg_cplex
2065 : type(atomdata_t) :: atom
2066 :
2067 : ! *************************************************************************
2068 : !write(ab_out,*)' prtdenmagsph : enter, rhomag(1,2)=',rhomag(1,2)
2069 :
2070 9802 : if(nspden==2)then
2071 1323 : rho_tot=rhomag(1,1) ; mag_coll=rhomag(1,2)
2072 1323 : if(cplex==2)then
2073 12 : rho_tot_im=rhomag(2,1) ; mag_coll_im=rhomag(2,2)
2074 : endif
2075 8479 : else if(nspden==4)then
2076 2156 : rho_tot=rhomag(1,1) ; mag_x=rhomag(1,2) ; mag_y=rhomag(1,3) ; mag_z=rhomag(1,4)
2077 2156 : if(cplex==2)then
2078 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)
2079 : endif
2080 : endif
2081 :
2082 9802 : if(option/=0)then
2083 :
2084 : !Printing
2085 9802 : sum_mag=zero
2086 9802 : sum_mag_x=zero
2087 9802 : sum_mag_y=zero
2088 9802 : sum_mag_z=zero
2089 9802 : sum_mag_x_im=zero
2090 9802 : sum_mag_y_im=zero
2091 9802 : sum_mag_z_im=zero
2092 9802 : sum_rho_up=zero
2093 9802 : sum_rho_dn=zero
2094 9802 : sum_rho_tot=zero
2095 :
2096 9802 : sum_mag_im=zero
2097 9802 : sum_mag_x_im=zero
2098 9802 : sum_mag_y_im=zero
2099 9802 : sum_mag_z_im=zero
2100 9802 : sum_rho_up_im=zero
2101 9802 : sum_rho_dn_im=zero
2102 9802 : sum_rho_tot_im=zero
2103 :
2104 : !Rotation matrices identity by default
2105 9802 : Rspin(:,:)=zero ; Rspin_t(:,:)=zero
2106 9802 : Rspin(1,1)=one ; Rspin(2,2)=one ; Rspin(3,3)=one
2107 9802 : Rspin_t(:,:)=Rspin(:,:)
2108 :
2109 : !Print spinaxis info only if axis not aligned with z
2110 9802 : if (abs(spinaxis(1))>tol8 .or. abs(spinaxis(2))>tol8) then
2111 1 : call geteuler(spinaxis, alpha, beta)
2112 1 : msg=' Spinaxis rotation information:'
2113 1 : write(msg, '(3a)' ) trim(msg),ch10,' ------------------------------'; call wrtout(units,msg)
2114 1 : write(msg, '(a,f12.6)') ' Alpha rotation angle around z-axis (degrees):', alpha * 180.0_dp / pi; call wrtout(units,msg)
2115 1 : write(msg, '(a,f12.6)') ' Beta rotation angle around y-axis (degrees): ', beta * 180.0_dp / pi; call wrtout(units,msg)
2116 1 : write(msg, '(a)') ' ---------------------------------------------------------'; call wrtout(units,msg)
2117 1 : call cart2spinaxis(alpha, beta, Rspin)
2118 1 : write(msg, '(a)') ' Rotation matrix from cartesian coordinate to spinaxis coordinate'; call wrtout(units,msg)
2119 1 : write(msg, '(3f14.8)') Rspin(1,1), Rspin(1,2), Rspin(1,3); call wrtout(units,msg)
2120 1 : write(msg, '(3f14.8)') Rspin(2,1), Rspin(2,2), Rspin(2,3); call wrtout(units,msg)
2121 1 : write(msg, '(3f14.8)') Rspin(3,1), Rspin(3,2), Rspin(3,3); call wrtout(units,msg)
2122 :
2123 13 : Rspin_t = transpose(Rspin)
2124 1 : write(msg, '(a)') ' Rotation matrix from spinaxis coordinate to cartesian coordinate'; call wrtout(units,msg)
2125 1 : write(msg, '(3f14.8)') Rspin_t(1,1), Rspin_t(1,2), Rspin_t(1,3); call wrtout(units,msg)
2126 1 : write(msg, '(3f14.8)') Rspin_t(2,1), Rspin_t(2,2), Rspin_t(2,3); call wrtout(units,msg)
2127 1 : write(msg, '(3f14.8)') Rspin_t(3,1), Rspin_t(3,2), Rspin_t(3,3); call wrtout(units,msg)
2128 1 : write(msg, '(a)') ' ----------------------------------------------------------------'
2129 1 : call wrtout(units,msg)
2130 : end if
2131 :
2132 9802 : mag_tot_cart=zero
2133 9802 : mag_tot_cart_im=zero
2134 9802 : if (nspden==4) then
2135 8624 : mag_spin = [mag_x,mag_y,mag_z]
2136 28028 : mag_tot_cart = matmul(Rspin_t,mag_spin)
2137 2156 : if (cplex==2) then
2138 44 : mag_spin = [mag_x_im,mag_y_im,mag_z_im]
2139 143 : mag_tot_cart_im = matmul(Rspin_t,mag_spin)
2140 : end if
2141 : end if
2142 :
2143 9802 : if(option==1 .or. option==11 .or. option==21) then
2144 :
2145 9802 : if(nspden==1) then
2146 6323 : if(option== 1)msg1=' Integrated electronic density in atomic spheres:'
2147 6323 : if(option==11)msg1=ch10//' Integrated potential residual in atomic spheres:'
2148 6323 : if(option==21)msg1=ch10//' Gradient with respect to target (=torque) :'
2149 6323 : write(msg, '(3a)' ) trim(msg1),ch10,' ------------------------------------------------'
2150 6323 : call wrtout(units,msg)
2151 6323 : if(ratsm>tol8)then
2152 221 : write(msg, '(a,f8.4,a)' ) ' Radius=ratsph(iatom), smearing ratsm=',ratsm,'.'
2153 221 : call wrtout(units,msg)
2154 : endif
2155 6323 : if(option== 1)then
2156 6213 : msg=' Atom Sphere_radius Integrated_density'
2157 6213 : if(present(ziontypat)) write(msg,'(a,a)')trim(msg),' Atomic charge'
2158 : endif
2159 6323 : if(option==11)msg=' Atom Sphere_radius Integrated_potresid'
2160 6323 : if(option==21)msg=' Atom Sphere_radius Torque'
2161 6323 : call wrtout(units,msg)
2162 25496 : do iatom=1,natom
2163 19173 : call atomdata_from_znucl(atom, znucl(typat(iatom)))
2164 19173 : write(msg, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),intgden(1,1,iatom)
2165 : !write(msg, '(i5,a3,f15.5,f20.8)' ) iatom,atom%symbol,ratsph(typat(iatom)),intgden(1,1,iatom)
2166 19173 : if(option==21)then
2167 : ! There is a change of sign to get the gradient wrt chrgat.
2168 10 : write(msg, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),-intgden(1,1,iatom)
2169 : !write(msg, '(i5,a3,f15.5,f20.8)' ) iatom,atom%symbol,ratsph(typat(iatom)),-intgden(1,1,iatom)
2170 : endif
2171 : !If option=1, print atomic charge
2172 19173 : if(option==1 .and. present(ziontypat))then
2173 10 : write(msg, '(a,f20.8)' ) trim(msg),ziontypat(typat(iatom))-intgden(1,1,iatom)
2174 : endif
2175 44669 : call wrtout(units,msg)
2176 : end do
2177 : endif
2178 :
2179 9802 : if(nspden==2 .or. nspden==4) then
2180 :
2181 3479 : if(option== 1)msg1=' Integrated electronic and magnetization densities in atomic spheres:'
2182 3479 : if(option==11)msg1=ch10//' Integrated potential residual in atomic spheres (scalar + magnetic field):'
2183 3479 : if(option==21)msg1=ch10//' Gradient with respect to target (=torque) (scalar + magnetic field):'
2184 3479 : write(msg, '(3a)' ) trim(msg1),ch10,' ------------------------------------------------------------------------------------------------------------'
2185 3479 : call wrtout(units,msg)
2186 :
2187 3479 : if(option== 1 .and. nspden==2) msg1='. Diff(up-dn)=approximate z local magnetic moment.'
2188 3479 : if(option== 1 .and. nspden==4) msg1='. mag(i)=approximate local magnetic moment.'
2189 3479 : if(option==11 .or. option==21) msg1='.'
2190 3479 : write(msg, '(a,f8.4,a)' ) ' Radius=ratsph(iatom), smearing ratsm=',ratsm,trim(msg1)
2191 3479 : call wrtout(units,msg)
2192 :
2193 3479 : if(option==1)then
2194 2231 : if(nspden==2) msg=' Atom Radius up_density dn_density Total(up+dn) Diff(up-dn)'
2195 2231 : if(nspden==4) msg=' Atom Radius Total density |mag| mag(x) mag(y) mag(z) mag(theta) mag(phi) '
2196 2231 : if(present(ziontypat))msg=trim(msg)//' Atomic charge'
2197 1248 : else if(option==11)then
2198 1199 : if(nspden==2) msg=' Atom Radius Potential B(z) up pot down pot'
2199 1199 : if(nspden==4) msg=' Atom Radius Potential |B| B(x) B(y) B(z) B(theta) B(phi) '
2200 49 : else if(option==21)then
2201 49 : if(nspden==2) msg=' Atom Radius grchrg T(z) up torque down torque'
2202 49 : if(nspden==4) msg=' Atom Radius Torque |T| T(x) T(y) T(z) T(theta) T(phi) '
2203 : endif
2204 3479 : if (cplex==2) then
2205 23 : msg_cplex=msg
2206 : else
2207 3456 : call wrtout(units,msg)
2208 : endif
2209 :
2210 : endif
2211 :
2212 9802 : if(nspden==2)then
2213 2658 : do icplex=1,cplex
2214 1335 : sum_mag=zero
2215 1335 : sum_rho_up=zero
2216 1335 : sum_rho_dn=zero
2217 1335 : sum_rho_tot=zero
2218 1335 : if (cplex==2) then
2219 24 : if (icplex==1) then
2220 12 : write(msg, '(a)') ' Real part of magnetization:'
2221 12 : call wrtout(units,msg)
2222 12 : call wrtout(units,msg_cplex)
2223 12 : elseif (icplex==2) then
2224 12 : write(msg, '(a)') '----------------------------'
2225 12 : call wrtout(units,msg)
2226 12 : write(msg, '(a)') ' Imaginary part of magnetization:'
2227 12 : call wrtout(units,msg)
2228 12 : call wrtout(units,msg_cplex)
2229 : endif
2230 : endif
2231 3839 : do iatom=1,natom
2232 2504 : call atomdata_from_znucl(atom, znucl(typat(iatom)))
2233 2504 : if(option/=21)then
2234 2480 : 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)
2235 : else
2236 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)
2237 : endif
2238 2504 : 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))
2239 2504 : if(option==1 .and. present(ziontypat)) then
2240 24 : write(msg, '(a,f14.6)') trim(msg),ziontypat(typat(iatom))-(intgden(icplex,1,iatom)+intgden(icplex,2,iatom))
2241 : end if
2242 2504 : call wrtout(units,msg)
2243 : ! Compute the sum of the magnetization
2244 2504 : sum_mag=sum_mag+intgden(icplex,1,iatom)-intgden(icplex,2,iatom)
2245 2504 : sum_rho_up=sum_rho_up+intgden(icplex,1,iatom)
2246 2504 : sum_rho_dn=sum_rho_dn+intgden(icplex,2,iatom)
2247 6343 : sum_rho_tot=sum_rho_tot+intgden(icplex,1,iatom)+intgden(icplex,2,iatom)
2248 : end do
2249 1335 : write(msg, '(a)') ' ------------------------------------------------------------------------------------------------------------'
2250 1335 : call wrtout(units,msg)
2251 1335 : write(msg, '(a,2f13.6,a,f12.6,a,f12.6)') ' Sum: ', sum_rho_up,sum_rho_dn,' ',sum_rho_tot,' ',sum_mag
2252 1335 : call wrtout(units,msg)
2253 :
2254 1335 : if(option==1)then
2255 1068 : write(msg, '(a,f14.6)') ' Total magnetization (from the atomic spheres): ', sum_mag
2256 1068 : call wrtout(units,msg)
2257 1068 : if (icplex==1) then
2258 1056 : write(msg, '(a,f14.6)') ' Total magnetization (exact up - dn): ', mag_coll
2259 12 : elseif (icplex==2) then
2260 12 : write(msg, '(a,f14.6)') ' Total magnetization (exact up - dn): ', mag_coll_im
2261 : endif
2262 1068 : call wrtout(units,msg)
2263 : endif
2264 1335 : write(msg, '(a)') ' '
2265 2658 : call wrtout(units,msg)
2266 : enddo
2267 :
2268 8479 : elseif(nspden==4) then
2269 :
2270 4323 : do icplex=1,cplex
2271 2167 : sum_mag_x=zero
2272 2167 : sum_mag_y=zero
2273 2167 : sum_mag_z=zero
2274 2167 : if (cplex==2) then
2275 22 : if (icplex==1) then
2276 11 : write(msg, '(a)') ' Real part of magnetization:'
2277 11 : call wrtout(units,msg)
2278 11 : call wrtout(units,msg_cplex)
2279 11 : elseif (icplex==2) then
2280 11 : write(msg, '(a)') '----------------------------'
2281 11 : call wrtout(units,msg)
2282 11 : write(msg, '(a)') ' Imaginary part of magnetization:'
2283 11 : call wrtout(units,msg)
2284 11 : call wrtout(units,msg_cplex)
2285 : endif
2286 : endif
2287 7661 : do iatom=1,natom
2288 71422 : mag_cart(:)=matmul(Rspin_t,intgden(icplex,2:4,iatom))
2289 5494 : call vcart2ylm(mag_cart, mag_r, mag_theta, mag_phi)
2290 5494 : call atomdata_from_znucl(atom, znucl(typat(iatom)))
2291 5494 : if(option/=21)then
2292 5402 : 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
2293 : else
2294 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
2295 : endif
2296 5494 : if(option==1 .and. present(ziontypat))&
2297 92 : & write(msg, '(a,f14.6)') trim(msg),ziontypat(typat(iatom))-intgden(icplex,1,iatom)
2298 5494 : call wrtout(units,msg)
2299 : ! Compute the sum of the magnetization in x, y and z directions
2300 5494 : sum_mag_x=sum_mag_x+mag_cart(1)
2301 5494 : sum_mag_y=sum_mag_y+mag_cart(2)
2302 13155 : sum_mag_z=sum_mag_z+mag_cart(3)
2303 : enddo
2304 8668 : vec = (/sum_mag_x,sum_mag_y,sum_mag_z /)
2305 2167 : call vcart2ylm(vec, sum_mag_r, sum_mag_theta, sum_mag_phi)
2306 2167 : write(msg, '(a)') ' ------------------------------------------------------------------------------------------------------------'
2307 2167 : call wrtout(units,msg)
2308 :
2309 6490 : if(option==1)then
2310 4724 : if (any(qgbt /= zero)) then
2311 8 : write(msg, '(a,f12.6,f12.6,f12.6)') ' Cell sum of sphere magnetization'
2312 8 : call wrtout(units,msg)
2313 : else
2314 1178 : write(msg, '(a)') ' Total magnetization (spheres) '
2315 1178 : call wrtout(units,msg)
2316 : endif
2317 1186 : write(msg, '(a,3f12.6)') ' (cart.coord.) ', sum_mag_x,sum_mag_y,sum_mag_z
2318 1186 : call wrtout(units,msg)
2319 1186 : write(msg, '(a,f12.6,a,2f12.6)') ' (sph.coord.) ', sum_mag_r," ",sum_mag_theta,sum_mag_phi
2320 1186 : call wrtout(units,msg)
2321 :
2322 4724 : if (any(qgbt /= zero)) then
2323 8 : write(msg, '(a,f10.6,f12.6,f12.6)') ' Integral of periodic magnetization'
2324 8 : call wrtout(units,msg)
2325 : else
2326 1178 : write(msg, '(a,f12.6,f12.6,4f12.6)') ' Total magnetization (exact) '
2327 1178 : call wrtout(units,msg)
2328 : endif
2329 1186 : if (icplex==1) then
2330 1175 : vec = mag_tot_cart
2331 1175 : call vcart2ylm(vec, exact_mag_r, exact_mag_theta, exact_mag_phi)
2332 1175 : write(msg, '(a,3f12.6)') ' (cart.coord.) ', mag_tot_cart(1),mag_tot_cart(2),mag_tot_cart(3)
2333 1175 : call wrtout(units,msg)
2334 1175 : write(msg, '(a,f12.6,a,2f12.6)') ' (sph.coord.) ', exact_mag_r," ",exact_mag_theta,exact_mag_phi
2335 1175 : call wrtout(units,msg)
2336 11 : elseif (icplex==2) then
2337 11 : vec = mag_tot_cart_im
2338 11 : call vcart2ylm(vec, exact_mag_r_im, exact_mag_theta_im, exact_mag_phi_im)
2339 11 : write(msg, '(a,3f12.6)') ' (cart.coord.) ', mag_tot_cart_im(1),mag_tot_cart_im(2),mag_tot_cart_im(3)
2340 11 : call wrtout(units,msg)
2341 11 : write(msg,'(a,f12.6,a,2f12.6)')' (sph.coord.) ', exact_mag_r_im," ", exact_mag_theta_im,exact_mag_phi_im
2342 11 : call wrtout(units,msg)
2343 : endif
2344 : endif
2345 : enddo
2346 : end if
2347 :
2348 : elseif(option==-1) then
2349 :
2350 0 : write(msg, '(2a)') ch10,' ---------------------------------------------------------------------------------------------------------------'
2351 0 : call wrtout(units,msg)
2352 :
2353 0 : if(nspden==1) then
2354 : write(msg, '(4a)' ) &
2355 0 : & ' Fermi level charge density n_f:',ch10,&
2356 0 : & ' ------------------------------------------------------------------------',ch10
2357 : else
2358 : write(msg, '(4a)' ) &
2359 0 : & ' Fermi level charge density n_f and magnetization m_f:',ch10,&
2360 0 : & ' ------------------------------------------------------------------------',ch10
2361 : end if
2362 0 : call wrtout(units,msg)
2363 :
2364 0 : if(cplex==1) then
2365 0 : write(msg, '(a,f13.8)') ' n_f = ',rho_tot
2366 : else
2367 0 : write(msg, '(a,f13.8,a,f13.8)') ' Re[n_f]= ', rho_tot," Im[n_f]= ",rho_tot_im
2368 : end if
2369 0 : call wrtout(units,msg)
2370 0 : if(nspden==2) then
2371 0 : if(cplex==1) then
2372 0 : write(msg, '(a,f13.8)') ' m_f = ', mag_coll
2373 : else
2374 0 : write(msg, '(a,f13.8,a,f13.8)') ' Re[m_f]= ', mag_coll," Im[m_f]= ",mag_coll_im
2375 : end if
2376 0 : call wrtout(units,msg)
2377 0 : elseif (nspden==4) then
2378 0 : write(msg, '(a,f13.8)') ' mx_f = ',mag_tot_cart(1)
2379 0 : call wrtout(units,msg)
2380 0 : write(msg, '(a,f13.8)') ' my_f = ',mag_tot_cart(2)
2381 0 : call wrtout(units,msg)
2382 0 : write(msg, '(a,f13.8)') ' mz_f = ',mag_tot_cart(3)
2383 0 : call wrtout(units,msg)
2384 : end if
2385 :
2386 0 : write(msg, '(3a)') ch10,' ------------------------------------------------------------------------',ch10
2387 0 : call wrtout(units,msg)
2388 :
2389 :
2390 : else if (option==2 .or. option==3 .or. option==4) then
2391 : ! Used in the DFPT case, option=idir+1
2392 :
2393 0 : if(abs(rho_tot)<1.0d-10) then
2394 0 : rho_tot=0
2395 : end if
2396 :
2397 0 : write(msg, '(2a)') ch10,' ---------------------------------------------------------------------------'
2398 0 : call wrtout(units,msg)
2399 :
2400 : ! if(nspden==1) then
2401 : ! write(msg, '(4a)' ) &
2402 : !& ' Integral of the first order density n^(1):',ch10,&
2403 : !& ' ------------------------------------------------------------------------',ch10
2404 : ! else
2405 : ! write(msg, '(4a)' ) &
2406 : !& ' Integrals of the first order density n^(1) and magnetization m^(1):',ch10,&
2407 : !& ' ------------------------------------------------------------------------',ch10
2408 : ! end if
2409 : ! call wrtout(units,msg)
2410 : !
2411 : ! if(cplex==1) then
2412 : ! write(msg, '(a,e16.8)') ' n^(1) = ', rho_tot
2413 : ! else
2414 : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[n^(1)] = ', rho_tot," Im[n^(1)] = ",rho_tot_im
2415 : ! end if
2416 : ! call wrtout(units,msg)
2417 : !
2418 : ! if(nspden==2) then
2419 : !
2420 : ! if(cplex==1) then
2421 : ! write(msg, '(a,e16.8)') ' m^(1) = ', mag_coll
2422 : ! else
2423 : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[m^(1)] = ', mag_coll," Im[m^(1)] = ",mag_coll_im
2424 : ! end if
2425 : ! call wrtout(units,msg)
2426 : !
2427 : ! elseif (nspden==4) then
2428 : ! if(cplex==1) then
2429 : ! write(msg, '(a,e16.8)') ' mx^(1) = ', mag_x
2430 : ! call wrtout(units,msg)
2431 : ! write(msg, '(a,e16.8)') ' my^(1) = ', mag_y
2432 : ! call wrtout(units,msg)
2433 : ! write(msg, '(a,e16.8)') ' mz^(1) = ', mag_z
2434 : ! call wrtout(units,msg)
2435 : ! else
2436 : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[mx^(1)]= ', mag_x, " Im[mx^(1)]= ", mag_x_im
2437 : ! call wrtout(units,msg)
2438 : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[my^(1)]= ', mag_y, " Im[my^(1)]= ", mag_y_im
2439 : ! call wrtout(units,msg)
2440 : ! write(msg, '(a,e16.8,a,e16.8)') ' Re[mz^(1)]= ', mag_z, " Im[mz^(1)]= ", mag_z_im
2441 : ! call wrtout(units,msg)
2442 : ! end if
2443 : ! end if
2444 : !
2445 : ! write(msg, '(3a)') ch10,' ------------------------------------------------------------------------',ch10
2446 : ! call wrtout(units,msg)
2447 :
2448 : end if
2449 :
2450 : end if ! option/=0
2451 :
2452 9802 : end subroutine prtdenmagsph
2453 : !!***
2454 :
2455 : !!****f* m_dens/radsmear
2456 : !! NAME
2457 : !! radsmear
2458 : !!
2459 : !! FUNCTION
2460 : !! As a function of the argument xarg (a positive number), return a function fsm that is zero
2461 : !! beyond some cut-off value xcut, one for xarg smaller than xcut-xsmear,
2462 : !! and interpolates smoothly between one and zero in the region from xcut-xsmear to xcut.
2463 : !! Also returns the derivative of this function, called dfsm.
2464 : !! The function fsm is twice differentiable at xcut (the first derivative is continuous, not the second),
2465 : !! and three times differentiable at xcut-xsmear (the second derivative is continuous, not the third).
2466 : !!
2467 : !!
2468 : !! INPUTS
2469 : !! xarg=argument of the function (should be positive)
2470 : !! xcut=largest value for which the function is non-zero
2471 : !! xsmear=defined the smearing region, between xcut-xsmear and xcut
2472 : !!
2473 : !! OUTPUT
2474 : !! fsm=value of the function
2475 : !! dfsm=derivative of the function with respect to xarg (zero, except in the smearing region).
2476 : !!
2477 : !! SOURCE
2478 :
2479 92223976 : subroutine radsmear(dfsm,fsm,xarg,xcut,xsmear)
2480 :
2481 : !Arguments ------------------------------------
2482 : !scalars
2483 : real(dp), intent(out) :: dfsm,fsm
2484 : real(dp), intent(in) :: xarg, xcut, xsmear
2485 :
2486 : !Local variables ------------------------------
2487 : !scalars
2488 : real(dp) :: xsmearinv,xx
2489 : !******************************************************************
2490 :
2491 92223976 : fsm = zero
2492 92223976 : dfsm=zero
2493 92223976 : if (xarg < xcut - xsmear - tol12) then
2494 90987178 : fsm = one
2495 1236798 : else if (xarg < xcut - tol12) then
2496 1236798 : xsmearinv=one/xsmear
2497 1236798 : xx = (xcut - xarg) * xsmearinv
2498 1236798 : fsm = xx**2*(3+xx*(1+xx*(-6+3*xx)))
2499 1236798 : dfsm = -(xx*(6+xx*(3+xx*(-24+15*xx))))*xsmearinv
2500 : end if
2501 :
2502 92223976 : end subroutine radsmear
2503 : !!***
2504 :
2505 : !!****f* ABINIT/printmagvtk
2506 : !! NAME
2507 : !! printmagvtk
2508 : !!
2509 : !! FUNCTION
2510 : !! Auxiliary routine for printing out magnetization density in VTK format.
2511 : !! Output file name is DEN.vtk
2512 : !!
2513 : !! INPUTS
2514 : !! mpi_enreg = information about adopted parallelization strategy
2515 : !! nspden = number of components of density matrix (possible values re 1,2, or 4)
2516 : !! nspden: 1 -> rho
2517 : !! nspden: 2 -> rho_up,rho_dwn
2518 : !! nspden: 4 -> rho,mx,my,mz
2519 : !! nfft = number of fft points per FFT processor
2520 : !! ngfft = full information about FFT mesh
2521 : !! rhor = density array stored in the memory of current FFT processor
2522 : !! rprimd = array of lattice vectors
2523 : !!
2524 : !! OUTPUT
2525 : !!
2526 : !! NOTES
2527 : !! At the moment this routine is mainly used for development and debugging
2528 : !! of gs and dfpt calculations with non-collinear spins. If needed, can be used
2529 : !! to print final density in vtk format.
2530 : !! IMPORTANT: implementation is thoroughly checked only for npspinor = 1,
2531 : !! for other case might need to change the part gathering
2532 : !! the FFT mesh info
2533 : !!
2534 : !! SOURCE
2535 :
2536 : subroutine printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,fname)
2537 :
2538 : !Arguments ------------------------------------
2539 : !scalars
2540 : type(MPI_type),intent(in) :: mpi_enreg
2541 : integer,intent(in) :: nfft,nspden,cplex
2542 : !arrays
2543 : integer,intent(in) :: ngfft(18)
2544 : real(dp),intent(in) :: rhor(cplex*nfft,nspden),rprimd(3,3)
2545 : character(len=*),intent(in) :: fname
2546 :
2547 : !Local variables-------------------------------
2548 : !scalars
2549 : integer :: denvtk,denxyz,denxyz_im,nfields
2550 : integer :: nx,ny,nz,nfft_tot
2551 : integer :: ii,jj,kk,ind,ispden
2552 : integer :: mpi_comm,mpi_head,mpi_rank,ierr
2553 : real(dp) :: rx,ry,rz
2554 : integer :: nproc_fft,ir
2555 : character(len=500) :: msg
2556 : character(len=10) :: outformat
2557 : character(len=fnlen) :: fname_vtk, fname_xyz, fname_xyz_re, fname_xyz_im
2558 : !arrays
2559 : real(dp),allocatable :: rhorfull(:,:)
2560 : ! *************************************************************************
2561 :
2562 : fname_vtk=adjustl(adjustr(fname)//".vtk")
2563 : fname_xyz=adjustl(adjustr(fname)//".xyz")
2564 : fname_xyz_re=adjustl(adjustr(fname)//"_re.xyz")
2565 : fname_xyz_im=adjustl(adjustr(fname)//"_im.xyz")
2566 : !write(std_out,*) ' Writing out .vtk file: ',fname_vtk
2567 : !write(std_out,*) ' Writing out .xyz file: ',fname_xyz
2568 :
2569 : !if 1 or two component density then write out either 1 or 2 scalar density fields
2570 : !if 4, then write one scalar field (density) and one vector field (magnetization density)
2571 : if(nspden/=4)then
2572 : nfields=nspden
2573 : else
2574 : nfields=2
2575 : end if
2576 :
2577 : nfields=nfields*cplex
2578 :
2579 : ! FFT mesh specifications: full grid
2580 : nx=ngfft(1) ! number of points along 1st lattice vector
2581 : ny=ngfft(2) ! number of points along 2nd lattice vector
2582 : nz=ngfft(3) ! number of points along 3rd lattice vector
2583 : 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)
2584 :
2585 :
2586 : ! Gather information about memory distribution
2587 : mpi_head=0
2588 : mpi_comm = mpi_enreg%comm_fft
2589 : mpi_rank = xmpi_comm_rank(mpi_comm)
2590 : nproc_fft=ngfft(10)
2591 :
2592 : ! Create array to host full FFT mesh
2593 : if(mpi_rank==mpi_head)then
2594 : ABI_MALLOC(rhorfull,(cplex*nfft_tot,nspden))
2595 : end if
2596 :
2597 : ! Fill in the full mesh
2598 : if(nproc_fft==1)then
2599 : rhorfull=rhor
2600 : else
2601 : do ir=1,nspden
2602 : call xmpi_gather(rhor(:,ir),cplex*nfft,rhorfull(:,ir),cplex*nfft,mpi_head,mpi_comm,ierr)
2603 : end do
2604 : end if
2605 :
2606 : if(mpi_rank==mpi_head)then
2607 :
2608 : ! Open the output vtk file
2609 : if (open_file(fname_vtk,msg,newunit=denvtk,status='replace',form='formatted') /=0) then
2610 : ABI_WARNING(msg)
2611 : RETURN
2612 : end if
2613 :
2614 : if(cplex==1) then
2615 : if (open_file(fname_xyz,msg,newunit=denxyz,status='replace',form='formatted') /=0) then
2616 : ABI_WARNING(msg)
2617 : RETURN
2618 : end if
2619 : else if (cplex==2) then
2620 : if (open_file(fname_xyz_re,msg,newunit=denxyz,status='replace',form='formatted') /=0) then
2621 : ABI_WARNING(msg)
2622 : RETURN
2623 : end if
2624 : if (open_file(fname_xyz_im,msg,newunit=denxyz_im,status='replace',form='formatted') /=0) then
2625 : ABI_WARNING(msg)
2626 : RETURN
2627 : end if
2628 : end if
2629 :
2630 : ! Write the header of the output vtk file
2631 : write(denvtk,"(a)") '# vtk DataFile Version 2.0'
2632 : write(denvtk,"(a)") 'Electron density components'
2633 : write(denvtk,"(a)") 'ASCII'
2634 : write(denvtk,"(a)") 'DATASET STRUCTURED_GRID'
2635 : write(denvtk,"(a,3i6)") 'DIMENSIONS ', nx,ny,nz
2636 : write(denvtk,"(a,i18,a)") 'POINTS ',nfft_tot,' double'
2637 :
2638 : if (nspden==1) then
2639 : outformat="(4e16.8)"
2640 : else if (nspden==2) then
2641 : outformat="(5e16.8)"
2642 : else
2643 : outformat="(7e16.8)"
2644 : end if
2645 :
2646 : ! Write out information about grid points
2647 : do kk=0,nz-1
2648 : do jj=0,ny-1
2649 : do ii=0,nx-1
2650 :
2651 : rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3)
2652 : ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3)
2653 : rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3)
2654 : write(denvtk,'(3f16.8)') rx,ry,rz !coordinates of the grid point
2655 : ind=1+ii+nx*(jj+ny*kk)
2656 : if (cplex==1) then
2657 : write(denxyz,outformat) rx,ry,rz,(rhorfull(ind,ispden),ispden=1,nspden)
2658 : else
2659 : write(denxyz,outformat) rx,ry,rz,(rhorfull(2*ind-1,ispden),ispden=1,nspden)
2660 : write(denxyz_im,outformat) rx,ry,rz,(rhorfull(2*ind ,ispden),ispden=1,nspden)
2661 : end if
2662 : end do
2663 : end do
2664 : end do
2665 :
2666 : if(cplex==1) then
2667 : close(denxyz)
2668 : else
2669 : close(denxyz)
2670 : close(denxyz_im)
2671 : end if
2672 :
2673 : ! Write out information about field defined on the FFT mesh
2674 : write(denvtk,"(a,i18)") 'POINT_DATA ',nfft_tot
2675 : write(denvtk,"(a,i6)") 'FIELD Densities ',nfields
2676 :
2677 :
2678 : ! Write out different fields depending on the number of density matrix components
2679 : if(nspden==1)then
2680 :
2681 : !single component, so just write out the density
2682 : if(cplex==1) then
2683 : write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
2684 : do kk=0,nz-1
2685 : do jj=0,ny-1
2686 : do ii=0,nx-1
2687 : ind=1+ii+nx*(jj+ny*kk)
2688 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2689 : end do
2690 : end do
2691 : end do
2692 : else
2693 : write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
2694 : do kk=0,nz-1
2695 : do jj=0,ny-1
2696 : do ii=0,nx-1
2697 : ind=2*(1+ii+nx*(jj+ny*kk))-1
2698 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2699 : end do
2700 : end do
2701 : end do
2702 : write(denvtk,"(a,i18,a)") 'Im_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))
2707 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2708 : end do
2709 : end do
2710 : end do
2711 : end if
2712 :
2713 : else if(nspden==2)then
2714 :
2715 : !two component, write the density for spin_up and spin_down channels
2716 : if(cplex==1) then
2717 :
2718 : write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
2719 : do kk=0,nz-1
2720 : do jj=0,ny-1
2721 : do ii=0,nx-1
2722 : ind=1+ii+nx*(jj+ny*kk)
2723 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2724 : end do
2725 : end do
2726 : end do
2727 : write(denvtk,"(a,i18,a)") 'mag 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)') 2*rhorfull(ind,2)-rhorfull(ind,1)
2733 : end do
2734 : end do
2735 : end do
2736 :
2737 : else
2738 :
2739 : write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
2740 : do kk=0,nz-1
2741 : do jj=0,ny-1
2742 : do ii=0,nx-1
2743 : ind=2*(1+ii+nx*(jj+ny*kk))-1
2744 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2745 : end do
2746 : end do
2747 : end do
2748 : write(denvtk,"(a,i18,a)") 'Im_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))
2753 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2754 : end do
2755 : end do
2756 : end do
2757 : write(denvtk,"(a,i18,a)") 'Re_mag 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))-1
2762 : write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
2763 : end do
2764 : end do
2765 : end do
2766 : write(denvtk,"(a,i18,a)") 'Im_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))
2771 : write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
2772 : end do
2773 : end do
2774 : end do
2775 :
2776 : end if
2777 :
2778 : else !here is the last option: nspden==4
2779 :
2780 : if(cplex==1) then
2781 :
2782 : !four component, write the density (scalar field) and magnetization density (vector field)
2783 : write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
2784 : do kk=0,nz-1
2785 : do jj=0,ny-1
2786 : do ii=0,nx-1
2787 : ind=1+ii+nx*(jj+ny*kk)
2788 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2789 : end do
2790 : end do
2791 : end do
2792 : write(denvtk,"(a,i18,a)") 'mag 3 ',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,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
2798 : end do
2799 : end do
2800 : end do
2801 :
2802 : else
2803 :
2804 : write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
2805 : do kk=0,nz-1
2806 : do jj=0,ny-1
2807 : do ii=0,nx-1
2808 : ind=2*(1+ii+nx*(jj+ny*kk))-1
2809 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2810 : end do
2811 : end do
2812 : end do
2813 : write(denvtk,"(a,i18,a)") 'Im_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))
2818 : write(denvtk,'(f16.8)') rhorfull(ind,1)
2819 : end do
2820 : end do
2821 : end do
2822 : write(denvtk,"(a,i18,a)") 'Re_mag 3 ',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))-1
2827 : write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
2828 : end do
2829 : end do
2830 : end do
2831 : write(denvtk,"(a,i18,a)") 'Im_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))
2836 : write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
2837 : end do
2838 : end do
2839 : end do
2840 :
2841 : end if
2842 :
2843 : end if ! nspden options condition
2844 :
2845 : close (denvtk)
2846 :
2847 : !clean up the gathered FFT mesh
2848 : ABI_FREE(rhorfull)
2849 :
2850 : end if
2851 :
2852 : end subroutine printmagvtk
2853 : !!***
2854 :
2855 : !!****f* ABINIT/calmaxdifmag
2856 : !! NAME
2857 : !! calmaxdifmag
2858 : !!
2859 : !! FUNCTION
2860 : !! Compute the maximum absolute value of the magnetization among all atoms,
2861 : !! and the maximum absolute change in magnetization between the current and previous SCF cycles.
2862 : !!
2863 : !! INPUTS
2864 : !! intgden(cplex,nspden,natom) = Integrated magnetic moments at the present SCF cycle.
2865 : !! intgden0(cplex,nspden,natom) = Integrated magnetic moments at the previous SCF cycle.
2866 : !! natom = Number of atoms.
2867 : !! nspden = Number of spin-density components (typically 1, 2, or 4).
2868 : !!
2869 : !! OUTPUT
2870 : !! maxmag = Maximum absolute value of magnetization among all atoms.
2871 : !! difmag = Maximum absolute change in magnetization between present and previous SCF cycle.
2872 : !!
2873 : !! SIDE EFFECTS
2874 : !! None
2875 : !!
2876 : !! NOTES
2877 : !! This routine can handle both real and complex spin densities.
2878 : !!
2879 : !! SOURCE
2880 :
2881 9102 : subroutine calmaxdifmag(cplex,intgden,intgden0,natom,nspden,maxmag,difmag)
2882 :
2883 : !Arguments ---------------------------------------------
2884 : integer,intent(in) :: natom,nspden,cplex
2885 : real(dp),intent(in) :: intgden(cplex,nspden,natom),intgden0(cplex,nspden,natom)
2886 : real(dp),intent(out)::maxmag,difmag
2887 : !Local variables ------------------------------
2888 : integer :: iatom
2889 : real(dp)::mag,mag0
2890 9102 : maxmag=zero;difmag=zero
2891 9102 : if (nspden==2 ) then
2892 19458 : do iatom=1,natom
2893 12992 : mag=intgden(1,1,iatom)-intgden(1,2,iatom)
2894 12992 : mag0=intgden0(1,1,iatom)-intgden0(1,2,iatom)
2895 12992 : maxmag=max(maxmag,abs(mag))
2896 12992 : difmag=max(difmag,abs(mag-mag0))
2897 19458 : if (cplex==2) then
2898 121 : mag=intgden(2,1,iatom)-intgden(2,2,iatom)
2899 121 : mag0=intgden(2,1,iatom)-intgden0(2,2,iatom)
2900 121 : maxmag=max(maxmag,abs(mag))
2901 121 : difmag=max(difmag,abs(mag-mag0))
2902 : endif
2903 : end do
2904 2636 : else if (nspden==4 ) then
2905 8392 : do iatom=1,natom
2906 23024 : mag =sqrt(sum(intgden(1,2:4,iatom)**2))
2907 23024 : mag0=sqrt(sum(intgden0(1,2:4,iatom)**2))
2908 5756 : maxmag=max(maxmag,abs(mag))
2909 5756 : difmag=max(difmag,abs(mag-mag0))
2910 8392 : if (cplex==2) then
2911 1240 : mag =sqrt(sum(intgden(2,2:4,iatom)**2))
2912 1240 : mag0=sqrt(sum(intgden0(2,2:4,iatom)**2))
2913 310 : maxmag=max(maxmag,abs(mag))
2914 310 : difmag=max(difmag,abs(mag-mag0))
2915 : endif
2916 : end do
2917 : endif
2918 9102 : if (maxmag < tol8) then
2919 2017 : maxmag=0
2920 2017 : if (difmag < tol8) difmag=0
2921 : endif
2922 9102 : end subroutine calmaxdifmag
2923 : !!***
2924 :
2925 : !!****f* m_dens/fatsph_recip
2926 : !! NAME
2927 : !! fatsph_recip
2928 : !!
2929 : !! FUNCTION
2930 : !! Compute the atomic spheres functions in reciprocal space as a product of
2931 : !! a Bessel function times a Gaussian smearing for the boundary. The functions
2932 : !! are subsequently Fourier transformed to real space.
2933 : !!
2934 : !! INPUTS
2935 : !! mpi_enreg=information about MPI parallelization
2936 : !! natom=number of atoms in cell.
2937 : !! nfft=(effective) number of FFT grid points (for this processor)
2938 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2939 : !! ntypat=number of atom types
2940 : !! ratsm=smearing width for ratsph
2941 : !! ratsph(ntypat)=radius of spheres around atoms
2942 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2943 : !! typat(natom)=type of each atom
2944 : !! xred(3,natom)=reduced dimensionless atomic coordinates
2945 : !!
2946 : !! OUTPUT
2947 : !! fatsph(nfft,natom)= functions defining the atomic spheres of integration in real space
2948 : !! fatsphi3(n1,n2,n3,natom)= Same functions with triple indexing
2949 : !!
2950 : !! SOURCE
2951 :
2952 0 : subroutine fatsph_recip(fatsph,fatsph3i,gmet,mpi_enreg,natom,nfft,ngfft,ntypat,&
2953 0 : & ratsm,ratsph,typat,ucvol,xred)
2954 :
2955 : !Arguments ------------------------------------
2956 : !scalars
2957 : integer,intent(in) :: natom,nfft,ntypat
2958 : real(dp),intent(in) :: ratsm,ucvol
2959 : type(MPI_type),intent(in) :: mpi_enreg
2960 : !arrays
2961 : integer,intent(in) :: ngfft(18),typat(natom)
2962 : real(dp),intent(in) :: gmet(3,3),ratsph(ntypat)
2963 : real(dp),intent(in) :: xred(3,natom)
2964 : real(dp),intent(out):: fatsph(nfft,natom)
2965 : real(dp),intent(out):: fatsph3i(ngfft(1),ngfft(2),ngfft(3),natom)
2966 :
2967 : !Local variables ------------------------------
2968 : !scalars
2969 : integer :: iatom
2970 : integer :: i1,i2,i3,id1,id2,id3,ig1,ig2,ig3,ii,ii1,n1,n2,n3,n4,n5,n6
2971 : real(dp) :: arg1,arg2,fac1,fac2,fac3,gq1,gq2,gq3,gcube
2972 : real(dp) :: gsquar,gmag,gmagrad,rad,sfr,sfi,widthsq
2973 : !arrays
2974 0 : integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
2975 0 : integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
2976 : real(dp) :: gq(3)
2977 0 : real(dp) :: work1(2,nfft),work2(nfft,1),work3(ngfft(1),ngfft(2),ngfft(3),1)
2978 :
2979 : !******************************************************************
2980 :
2981 : !Geometric parameters
2982 0 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
2983 0 : n4=ngfft(4);n5=ngfft(5);n6=ngfft(6)
2984 0 : id1=n1/2+2
2985 0 : id2=n2/2+2
2986 0 : id3=n3/2+2
2987 0 : widthsq=ratsm**2
2988 :
2989 : !Get the distrib associated with this fft_grid
2990 0 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2991 :
2992 0 : do iatom=1, natom
2993 :
2994 0 : ii=0
2995 0 : work1(:,:)=zero
2996 : !G=0 term
2997 0 : rad=ratsph(typat(iatom))
2998 0 : work1(1,1)=four_pi*rad**3/(three*ucvol)
2999 0 : do i3=1,n3
3000 0 : ig3=i3-(i3/id3)*n3-1
3001 0 : gq3=dble(ig3)
3002 0 : gq(3)=gq3
3003 0 : do i2=1,n2
3004 0 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
3005 0 : ig2=i2-(i2/id2)*n2-1
3006 0 : gq2=dble(ig2)
3007 0 : gq(2)=gq2
3008 :
3009 : ! Note the lower limit of the next loop
3010 0 : ii1=1
3011 0 : if(i3==1 .and. i2==1 .and. ig2==0 .and. ig3==0)then
3012 0 : ii1=2
3013 0 : ii=ii+1
3014 : end if
3015 0 : do i1=ii1,n1
3016 0 : ig1=i1-(i1/id1)*n1-1
3017 0 : gq1=dble(ig1)
3018 0 : gq(1)=gq1
3019 0 : ii=ii+1
3020 :
3021 0 : gsquar=gsq_vl3(gq1,gq2,gq3)
3022 0 : gmag=sqrt(gsquar)
3023 0 : gcube=gmag*gsquar
3024 0 : gmagrad=two_pi*gmag*rad
3025 0 : arg1=-gsquar*pi**2*widthsq
3026 0 : arg2=two_pi*dot_product(xred(:,iatom),gq)
3027 :
3028 0 : fac1=two/(two_pi**2*gcube*ucvol)
3029 0 : fac2=sin(gmagrad)-gmagrad*cos(gmagrad)
3030 0 : fac3=exp(arg1)
3031 0 : sfr=cos(arg2)
3032 0 : sfi=-sin(arg2)
3033 :
3034 0 : work1(1,ii)=fac1*fac2*fac3*sfr
3035 0 : work1(2,ii)=fac1*fac2*fac3*sfi
3036 :
3037 : end do
3038 : end if
3039 : end do
3040 : end do
3041 :
3042 : ! Transform to real space
3043 0 : call fourdp(1,work1,work2,1,mpi_enreg,nfft,1,ngfft,0)
3044 0 : fatsph(:,iatom)=work2(:,1)
3045 :
3046 0 : call fftpac(1,mpi_enreg,1,n1,n2,n3,n1,n2,n3,ngfft,work2,work3,2)
3047 0 : fatsph3i(:,:,:,iatom)=work3(:,:,:,1)
3048 :
3049 : end do !iatom
3050 :
3051 : contains
3052 :
3053 0 : function gsq_vl3(g1,g2,g3)
3054 :
3055 : real(dp) :: gsq_vl3
3056 : real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions
3057 : !Define G^2 based on G space metric gmet.
3058 : gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+&
3059 : & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+&
3060 0 : & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1)
3061 0 : end function gsq_vl3
3062 :
3063 : end subroutine fatsph_recip
3064 : !!***
3065 :
3066 : !!****f* m_dens/magmom_to_d2
3067 : !! NAME
3068 : !! magmom_to_d2
3069 : !!
3070 : !! FUNCTION
3071 : !! Incorporates the magnetic moments in the ddb files as second
3072 : !! order energy derivatives with respect to (ipert,idir) and a
3073 : !! Zeeman field. Both total, i.e., response to a uniform Zeeman
3074 : !! field (ipert=natom+5) and local (ipert=natom+11+1:2*natom+11)
3075 : !! magnetic moments are considered. The Zeeman field directions
3076 : !! are passed in Cartesian format.
3077 : !! It also incorporates the second-order energy derivatives involving
3078 : !! a scalar-potential perturbation from the charge-induced by
3079 : !! another (ipert,idir). TODO: Check the consistency of signs in this case.
3080 : !!
3081 : !! INPUTS
3082 : !! blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
3083 : !! idir=direction of the perturbation
3084 : !! intgden(cplex,nspden, natom)=integrated rhor or potential residual, for each atom in a sphere of radius ratsph.
3085 : !! Representation differs according to nspden :
3086 : !! if nspden=1, total density
3087 : !! if nspden=2, spin up, then spin down
3088 : !! if nspden=4, total density, then mag_x, mag_y, mag_z
3089 : !! ipert=type of perturbation
3090 : !! mpert=maximum number of perturbations
3091 : !! natom=number of atoms in cell.
3092 : !! nspden=number of spin-density components
3093 : !! rhomag(2,nspden)=integral of charge or magnetization over the whole cell (also taking into account a possible imaginary part for DFPT).
3094 : !!
3095 : !! OUTPUT
3096 : !! d2lo(2,3,mpert,3,mpert)= Local contributions to the second-order energy functional.
3097 : !!
3098 : !! SOURCE
3099 :
3100 12 : subroutine magmom_to_d2(blkflg,cplex,d2lo,idir,intgden,ipert,mpert,natom,nspden,rhomag)
3101 :
3102 : !Arguments ---------------------------------------------
3103 : !scalars
3104 : integer,intent(in) :: cplex,idir,ipert,mpert,natom,nspden
3105 : !arrays
3106 : integer,intent(inout) :: blkflg(3,mpert,3,mpert)
3107 : real(dp),intent(in) :: intgden(cplex,nspden,natom)
3108 : real(dp),intent(in) :: rhomag(2,nspden)
3109 : real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert)
3110 : !Local variables ------------------------------
3111 : !scalars
3112 : integer :: iatom
3113 :
3114 : ! *************************************************************************
3115 :
3116 : ! We store in DDB the second-order energy derivatives, hence the negative
3117 : ! sign applied to the induced magnetic moments.
3118 :
3119 : ! Incorporate total charge and magnetic moments
3120 12 : if (nspden==2) then
3121 0 : blkflg(1,natom+6,idir,ipert)= 1
3122 0 : d2lo(1,1,natom+6,idir,ipert)= -rhomag(1,1)
3123 0 : if (cplex==2) d2lo(2,1,natom+6,idir,ipert)= rhomag(2,1)
3124 0 : blkflg(3,natom+5,idir,ipert)= 1
3125 0 : d2lo(1,3,natom+5,idir,ipert)= -half*rhomag(1,2)
3126 0 : if (cplex==2) d2lo(2,3,natom+5,idir,ipert)= -half*rhomag(2,2)
3127 12 : else if (nspden==4) then
3128 12 : blkflg(1,natom+6,idir,ipert)=1
3129 12 : d2lo(1,1,natom+6,idir,ipert)= -rhomag(1,1)
3130 12 : if (cplex==2) d2lo(2,1,natom+6,idir,ipert)= rhomag(2,1)
3131 48 : blkflg(1:3,natom+5,idir,ipert)=1
3132 48 : d2lo(1,1:3,natom+5,idir,ipert)= -half*rhomag(1,2:4)
3133 36 : if (cplex==2) d2lo(2,1:3,natom+5,idir,ipert)= -half*rhomag(2,2:4)
3134 : end if
3135 :
3136 : ! Incorporate local magnetic moments
3137 12 : if (nspden==2) then
3138 0 : do iatom= 1, natom
3139 0 : blkflg(3,natom+11+iatom,idir,ipert)= 1
3140 0 : d2lo(1,3,natom+11+iatom,idir,ipert)= -half*intgden(1,2,iatom)
3141 0 : if (cplex==2) d2lo(2,3,natom+11+iatom,idir,ipert)= -half*intgden(2,2,iatom)
3142 : end do
3143 12 : else if (nspden==4) then
3144 84 : do iatom= 1, natom
3145 288 : blkflg(1:3,natom+11+iatom,idir,ipert)= 1
3146 288 : d2lo(1,1:3,natom+11+iatom,idir,ipert)= -half*intgden(1,2:4,iatom)
3147 228 : if (cplex==2) d2lo(2,1:3,natom+11+iatom,idir,ipert)= -half*intgden(2,2:4,iatom)
3148 : end do
3149 : end if
3150 :
3151 12 : end subroutine magmom_to_d2
3152 : !!***
3153 :
3154 0 : end module m_dens
3155 : !!***
|
