Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_types_jij
8 :
9 : USE m_types
10 : USE m_types_forcetheo
11 : USE m_judft
12 : #ifdef CPP_MPI
13 : USE mpi
14 : #endif
15 : TYPE,EXTENDS(t_forcetheo) :: t_forcetheo_jij
16 : INTEGER :: loopindex,no_loops
17 : INTEGER,ALLOCATABLE :: q_index(:),iatom(:),jatom(:)
18 : LOGICAL,ALLOCATABLE :: phase2(:)
19 :
20 : REAL,ALLOCATABLE:: qvec(:,:)
21 : REAL :: thetaj
22 : REAL,ALLOCATABLE:: evsum(:)
23 : CONTAINS
24 : PROCEDURE :: start =>jij_start
25 : PROCEDURE :: next_job=>jij_next_job
26 : PROCEDURE :: eval =>jij_eval
27 : PROCEDURE :: postprocess => jij_postprocess
28 : PROCEDURE :: init => jij_init !not overloaded
29 : PROCEDURE :: dist => jij_dist !not overloaded
30 : END TYPE t_forcetheo_jij
31 :
32 : CONTAINS
33 :
34 :
35 :
36 :
37 0 : SUBROUTINE jij_init(this,qvec,thetaj,atoms)
38 : USE m_types_setup
39 : USE m_constants
40 : IMPLICIT NONE
41 : CLASS(t_forcetheo_jij),INTENT(INOUT):: this
42 : REAL,INTENT(in) :: qvec(:,:),thetaj
43 : TYPE(t_atoms),INTENT(IN) :: atoms
44 :
45 : INTEGER:: n,na,ni,nj,j
46 : REAL,PARAMETER:: eps=1E-5
47 :
48 0 : this%l_needs_vectors=.false.
49 :
50 0 : this%qvec=qvec
51 0 : this%thetaj=thetaj
52 :
53 : !Max no of loops...
54 0 : n=atoms%nat**2*SIZE(this%qvec,2)+1
55 0 : ALLOCATE(this%q_index(n),this%iatom(n),this%jatom(n),this%phase2(n))
56 :
57 :
58 : !now construct the loops
59 0 : this%no_loops=0
60 0 : DO n=1,SIZE(this%qvec,2)
61 0 : DO ni=1,atoms%ntype
62 0 : IF (ABS(atoms%bmu(ni))<eps.and..not.atoms%econf(ni)%is_polarized()) CYCLE !no magnetic atom
63 0 : DO nj=ni,atoms%ntype
64 0 : IF (ABS(atoms%bmu(nj))<eps.and..not.atoms%econf(nj)%is_polarized()) CYCLE !no magnetic atom
65 0 : DO j=1,MERGE(1,2,ni==nj) !phase factor
66 : !new config found
67 0 : this%no_loops=this%no_loops+1
68 0 : this%q_index(this%no_loops)=n
69 0 : this%iatom(this%no_loops)=ni
70 0 : this%jatom(this%no_loops)=nj
71 0 : this%phase2(this%no_loops)=(j==2)
72 : ENDDO
73 : END DO
74 : END DO
75 : END DO
76 0 : if (this%no_loops == 0) then
77 0 : call judft_error("No configurations for Jij calculation set. Make sure you have magnetic atoms in the unit cell")
78 : endif
79 :
80 0 : ALLOCATE(this%evsum(this%no_loops))
81 0 : this%evsum=0
82 0 : END SUBROUTINE jij_init
83 :
84 :
85 0 : SUBROUTINE jij_dist(this,fmpi)
86 : USE m_types_mpi
87 : IMPLICIT NONE
88 : CLASS(t_forcetheo_jij),INTENT(INOUT):: this
89 : TYPE(t_mpi),INTENT(in):: fmpi
90 :
91 : INTEGER:: i,q,ierr
92 : #ifdef CPP_MPI
93 0 : call judft_error("jij has to be parallelized")
94 : #endif
95 0 : END SUBROUTINE jij_dist
96 :
97 :
98 :
99 0 : SUBROUTINE jij_start(this,potden,l_io)
100 : USE m_types_potden
101 : IMPLICIT NONE
102 : CLASS(t_forcetheo_jij),INTENT(INOUT):: this
103 : TYPE(t_potden) ,INTENT(INOUT) :: potden
104 : LOGICAL,INTENT(IN) :: l_io
105 0 : this%loopindex=0
106 0 : CALL this%t_forcetheo%start(potden,l_io) !call routine of basis type
107 0 : END SUBROUTINE jij_start
108 :
109 0 : LOGICAL FUNCTION jij_next_job(this,fmpi,lastiter,atoms,noco,nococonv)
110 : USE m_types_setup
111 : USE m_xmlOutput
112 : USE m_constants
113 : USE m_types_nococonv
114 : USE m_types_mpi
115 : IMPLICIT NONE
116 : CLASS(t_forcetheo_jij),INTENT(INOUT):: this
117 : TYPE(t_mpi), INTENT(IN) :: fmpi
118 : LOGICAL,INTENT(IN) :: lastiter
119 : TYPE(t_atoms),INTENT(IN) :: atoms
120 : TYPE(t_noco),INTENT(IN) :: noco
121 : !Stuff that might be modified...
122 : TYPE(t_nococonv),INTENT(INOUT) :: nococonv
123 :
124 : !locals
125 : INTEGER:: n
126 : CHARACTER(LEN=12):: attributes(2)
127 :
128 0 : IF (.NOT.lastiter) THEN
129 0 : jij_next_job=this%t_forcetheo%next_job(fmpi,lastiter,atoms,noco,nococonv)
130 0 : RETURN
131 : ENDIF
132 :
133 : !OK, now we start the JIJ-loop
134 0 : this%l_in_forcetheo_loop = .true.
135 0 : this%loopindex=this%loopindex+1
136 0 : jij_next_job=(this%loopindex<=this%no_loops) !still loops to do...
137 0 : IF (.NOT.jij_next_job) RETURN
138 :
139 : ! Now set the noco-variable accordingly...
140 :
141 0 : nococonv%qss=this%qvec(:,this%q_index(this%loopindex))
142 :
143 : !c Determines the cone angles and phase shifts of the spin
144 : !c vectors on magnetic atoms for the calculation of the
145 : !c interaction constants Jij from the Heisenberg model
146 : !c M. Lezaic 04
147 :
148 0 : nococonv%alph=0.0;nococonv%beta=0.0
149 0 : nococonv%beta(this%iatom(this%loopindex))=this%thetaj
150 0 : IF (this%phase2(this%loopindex)) nococonv%alph(this%iatom(this%loopindex))=pi_const*0.5
151 0 : nococonv%beta(this%jatom(this%loopindex))=this%thetaj
152 :
153 : !rotate according to q-vector
154 0 : DO n = 1,atoms%ntype
155 0 : nococonv%alph(n) = nococonv%alph(n) + tpi_const*DOT_PRODUCT(nococonv%qss,atoms%taual(:,atoms%firstAtom(n)))
156 : ENDDO
157 :
158 0 : IF (.NOT.this%l_io) RETURN
159 0 : IF (fmpi%irank .EQ. 0) THEN
160 0 : IF (this%loopindex.NE.1) CALL closeXMLElement('Forcetheorem_Loop')
161 0 : WRITE(attributes(1),'(a)') 'JIJ'
162 0 : WRITE(attributes(2),'(i5)') this%loopindex
163 0 : CALL openXMLElementPoly('Forcetheorem_Loop',(/'calculationType','No '/),attributes)
164 : END IF
165 : END FUNCTION jij_next_job
166 :
167 0 : SUBROUTINE jij_postprocess(this)
168 : USE m_xmlOutput
169 : IMPLICIT NONE
170 : CLASS(t_forcetheo_jij),INTENT(INOUT):: this
171 :
172 : !Locals
173 : INTEGER:: n
174 : CHARACTER(LEN=18):: attributes(6)
175 :
176 0 : IF (this%loopindex==0) RETURN
177 :
178 0 : IF (.NOT.this%l_io) RETURN
179 :
180 : !Now output the results
181 0 : call closeXMLElement('Forcetheorem_Loop')
182 0 : attributes = ''
183 0 : WRITE(attributes(1),'(i5)') this%no_loops
184 0 : WRITE(attributes(2),'(a)') 'Htr'
185 0 : CALL openXMLElement('Forcetheorem_JIJ',(/'Configs','units '/),attributes(:2))
186 0 : DO n=1,this%no_loops
187 0 : WRITE(attributes(1),'(i5)') n
188 0 : WRITE(attributes(2),'(3(f5.3,1x))') this%qvec(:,this%q_index(n))
189 0 : WRITE(attributes(3),'(i0)') this%iatom(n)
190 0 : WRITE(attributes(4),'(i0)') this%jatom(n)
191 0 : WRITE(attributes(5),'(l1)') this%phase2(n)
192 0 : WRITE(attributes(6),'(f15.8)') this%evsum(n)
193 : CALL writeXMLElementForm('Config',(/'n ','q ','iatom ','jatom ','phase ','ev-sum'/),attributes,&
194 0 : RESHAPE((/1,1,5,5,5,6,4,18,4,4,2,15/),(/6,2/)))
195 : ENDDO
196 0 : CALL closeXMLElement('Forcetheorem_JIJ')
197 0 : CALL judft_end("Forcetheorem: Jij")
198 : END SUBROUTINE jij_postprocess
199 :
200 :
201 0 : FUNCTION jij_eval(this,eig_id,atoms,kpts,sym,&
202 : cell,noco,nococonv, input,fmpi, enpara,v,results)RESULT(skip)
203 : USE m_types
204 : USE m_ssomat
205 : IMPLICIT NONE
206 : CLASS(t_forcetheo_jij),INTENT(INOUT):: this
207 : LOGICAL :: skip
208 : !Stuff that might be used...
209 : TYPE(t_mpi),INTENT(IN) :: fmpi
210 :
211 :
212 : TYPE(t_input),INTENT(IN) :: input
213 : TYPE(t_noco),INTENT(IN) :: noco
214 : TYPE(t_nococonv),INTENT(IN) :: nococonv
215 : TYPE(t_sym),INTENT(IN) :: sym
216 : TYPE(t_cell),INTENT(IN) :: cell
217 : TYPE(t_kpts),INTENT(IN) :: kpts
218 : TYPE(t_atoms),INTENT(IN) :: atoms
219 : TYPE(t_enpara),INTENT(IN) :: enpara
220 : TYPE(t_potden),INTENT(IN) :: v
221 : TYPE(t_results),INTENT(IN) :: results
222 : INTEGER,INTENT(IN) :: eig_id
223 0 : skip=.FALSE.
224 0 : IF (this%loopindex==0) RETURN
225 :
226 0 : this%evsum(this%loopindex)=results%seigv
227 0 : skip=.TRUE.
228 0 : END FUNCTION jij_eval
229 :
230 :
231 :
232 0 : SUBROUTINE priv_analyse_data()
233 : !-------------------------------------------------------------------
234 : ! calculates the
235 : ! coupling constants J for all the couples of magnetic
236 : ! atoms, for a given number of neighbouring shells
237 : ! M. Lezaic 04
238 : !-------------------------------------------------------------------
239 :
240 : USE m_constants
241 0 : PRINT *,"jcoef2 has still to be reimplemented"
242 : #ifdef CPP_NEVER
243 : USE m_nshell
244 : IMPLICIT NONE
245 :
246 : c .. Scalar arguments ..
247 :
248 : INTEGER, INTENT (IN) :: ntypd,nmagn,nqpt,mtypes
249 : INTEGER, INTENT (IN) :: natd,nop
250 : INTEGER, INTENT (INOUT) :: nsh
251 : REAL, INTENT (IN) :: thetaJ
252 : LOGICAL, INTENT (IN) :: invs,film
253 :
254 : c .. Array arguments ..
255 :
256 : INTEGER, INTENT (IN) :: mrot(3,3,nop)
257 : INTEGER, INTENT (IN) :: nmagtype(ntypd),magtype(ntypd),neq(ntypd)
258 : REAL, INTENT (IN) :: taual(3,natd)
259 : REAL, INTENT (IN) :: amat(3,3)
260 :
261 : c .. Temporary..
262 : INTEGER :: nshort !The number of non-zero interactions for a calculation
263 : c with the least square method
264 : REAL Tc
265 : REAL, ALLOCATABLE :: Cmat(:,:),DelE(:),work(:)
266 : INTEGER lwork
267 :
268 : c .. Local scalars ..
269 :
270 : INTEGER n,nn,nnn,mu,nu,iatom,nqvect,atsh,phi,i,ii,wrJ
271 : INTEGER qcount,info,imt,limit,remt,sneq,nqcomp
272 : REAL qn,sqsin,alph,beta,tpi,scp,J,rseig
273 : REAL enmin,enmax,zcoord,deltaz
274 : ! REAL wei !(old version)
275 : INTEGER, PARAMETER :: nmax=10,dims=(2*nmax+1)**3,shmax=192
276 : REAL, PARAMETER :: tol=0.000001
277 :
278 : c ..Local Arrays..
279 :
280 : INTEGER w(nop,nqpt-1),itype(nmagn)
281 : INTEGER nat(dims),ierr(3)
282 : REAL seigv(nmagn,nmagn,nqpt-1,2),q(3,nop,nqpt-1),qss(3),Jw(shmax)
283 : REAL tauC1(3),tauC2(3),seigv0(nmagn,nmagn,2)
284 : REAL ReJq(nmagn,nmagn,nqpt-1),ImJq(nmagn,nmagn,nqpt-1)
285 : REAL M(nmagn),qssmin(3),qssmax(3),t(3),R(3,shmax,dims),lenR(dims)
286 : REAL Dabsq(3),divi(3)
287 : INTEGER IDabsq(3)
288 :
289 : c .. Intrinsic Functions ..
290 :
291 : INTRINSIC cos,sin
292 :
293 : c .. External Subroutines ..
294 :
295 : EXTERNAL CPP_LAPACK_dgels
296 :
297 : c-------------------------------------------------------------------
298 : OPEN(116,file='qptsinfo',status='unknown')
299 : OPEN(115,file='jconst',status='unknown')
300 : w=0
301 : nqvect=0
302 : nshort=0
303 : ReJq=0.0
304 : ImJq=0.0
305 : sqsin=(sin(thetaJ))**2
306 : tpi = 2.0 * pimach()
307 : limit=nmagn-1
308 : IF (nmagn.gt.mtypes) limit=mtypes
309 :
310 : IF(nsh.LT.0)THEN
311 : c... Setup for a calculation using least squares method
312 : WRITE(oUnit,*) 'Jij calculated with the least squares method'
313 : nsh=-nsh
314 : nshort=nsh
315 : ENDIF
316 :
317 : DO n=1,nqpt
318 : qcount=n-1
319 : mu=1
320 : DO imt=1,mtypes
321 : DO nu=mu,nmagn
322 : DO phi=1,2
323 : READ(114,*)
324 : READ(114,5000) itype(mu),alph,beta,M(mu)
325 : READ(114,5000) itype(nu),alph,beta,M(nu)
326 : ! READ(114,5001) qss(1),qss(2),qss(3),wei, rseig !(old version)
327 : READ(114,5001) qss(1),qss(2),qss(3),rseig
328 : 5000 FORMAT(3x,i4,4x,f14.10,1x,f14.10,4x,f8.5)
329 : !5001 FORMAT(4(f14.10,1x),f20.10) !(old version)
330 : 5001 FORMAT(3(f14.10,1x),f20.10)
331 :
332 : c... Aquire information on the maximal and minimal calculated energy
333 : IF((n*mu*nu*phi).EQ.1)THEN
334 : qssmin=qss
335 : enmin=rseig
336 : qssmax=qss
337 : enmax=rseig
338 : ELSE
339 : IF(enmin.GE.rseig)THEN
340 : enmin=rseig
341 : qssmin=qss
342 : ENDIF
343 : IF(enmax.LE.rseig)THEN
344 : enmax=rseig
345 : qssmax=qss
346 : ENDIF
347 : ENDIF
348 :
349 : IF(n.EQ.1)THEN
350 : seigv0(mu,nu,phi)=rseig ! reference:energy of the collinear state
351 : ELSE
352 : seigv(mu,nu,qcount,phi)=rseig ! energies of spin-spirals
353 : ENDIF
354 : qn = ( qss(1)**2 + qss(2)**2 + qss(3)**2 )
355 : IF((mu.EQ.nu).OR.(invs.AND.(qn.GT.tol)))
356 : & GOTO 33
357 : ENDDO !phi
358 : 33 CONTINUE
359 : ENDDO !nu
360 : mu=mu+nmagtype(imt)
361 : ENDDO !imt
362 : IF(n.EQ.1)THEN
363 : IF(qn.GT.tol) THEN
364 : WRITE(*,*) 'jcoff2: The first energy in jenerg file should correspond
365 : &to qss=0!'
366 : #ifdef CPP_MPI
367 : CALL MPI_ABORT(MPI_COMM_WORLD,1,ierr)
368 : #endif
369 : STOP
370 : ENDIF
371 : ELSE
372 : WRITE(116,*) qcount
373 : c...
374 : c... Apply all the symmetry operations to the q-vectors from the irreducible wedge
375 : c...
376 : DO nn=1,nop
377 : q(1,nn,qcount)=qss(1)*mrot(1,1,nn) + qss(2)*mrot(2,1,nn) +
378 : + qss(3)*mrot(3,1,nn)
379 : q(2,nn,qcount)=qss(1)*mrot(1,2,nn) + qss(2)*mrot(2,2,nn) +
380 : + qss(3)*mrot(3,2,nn)
381 : q(3,nn,qcount)=qss(1)*mrot(1,3,nn) + qss(2)*mrot(2,3,nn) +
382 : + qss(3)*mrot(3,3,nn)
383 : w(nn,qcount)=1
384 : c...
385 : c... Eliminate the q-vectors which are repeating and for each q remove the -q
386 : c...
387 : DO nnn=1,nn-1
388 : DO ii=1,2
389 : Dabsq(:)=ABS(q(:,nn,qcount)+((-1)**ii)*q(:,nnn,qcount))
390 : IDabsq(:)=NINT(Dabsq(:))
391 : divi(:)=ABS(Dabsq(:)/FLOAT(IDabsq(:))-1.0)
392 : IF(((Dabsq(1).LT.tol).OR.(divi(1).LT.tol)).AND.
393 : & ((Dabsq(2).LT.tol).OR.(divi(2).LT.tol)).AND.
394 : & ((Dabsq(3).LT.tol).OR.(divi(3).LT.tol)))THEN
395 : w(nn,qcount)=0
396 : GOTO 44
397 : ENDIF
398 : ENDDO ! ii
399 : ENDDO !nnn
400 : nqvect=nqvect+1
401 : WRITE(116,5502) mrot(1,1,nn),mrot(1,2,nn),mrot(1,3,nn)
402 : WRITE(116,5502) mrot(2,1,nn),mrot(2,2,nn),mrot(2,3,nn)
403 : WRITE(116,5502) mrot(3,1,nn),mrot(3,2,nn),mrot(3,3,nn)
404 : WRITE(116,5002) nqvect,q(1,nn,qcount),q(2,nn,qcount),
405 : & q(3,nn,qcount)
406 : 5002 FORMAT(i6,3(f14.10,1x))
407 : 5502 FORMAT(3(i3))
408 : 44 CONTINUE
409 : ENDDO !nn
410 : c... Now calculate Jq=Re(Jq)+i*Im(Jq)
411 : mu=1
412 : DO imt=1,mtypes
413 : ReJq(mu,mu,qcount)=-2.0*(seigv(mu,mu,qcount,1)
414 : & -seigv0(mu,mu,1))/(M(mu)*M(mu)*sqsin)
415 : ImJq(mu,mu,qcount)=0.0
416 : DO remt=mu+1,mu+nmagtype(imt)-1
417 : ReJq(remt,remt,qcount)=ReJq(mu,mu,qcount)
418 : ImJq(remt,remt,qcount)=0.0
419 : ENDDO!remt
420 : mu=mu+nmagtype(imt)
421 : ENDDO !imt
422 : mu=1
423 : DO imt=1,limit
424 : DO nu=mu+1,nmagn
425 : ReJq(mu,nu,qcount)=((seigv0(mu,nu,2)-
426 : & seigv(mu,nu,qcount,1))/(M(mu)*M(nu)*sqsin))
427 : & -(0.5*M(mu)*ReJq(mu,mu,qcount)/M(nu))-
428 : & (0.5*M(nu)*ReJq(nu,nu,qcount)/M(mu))
429 : IF(invs)THEN
430 : ImJq(mu,nu,qcount)=0.0
431 : ELSE
432 : ImJq(mu,nu,qcount)=((seigv(mu,nu,qcount,2)
433 : & -seigv(mu,nu,qcount,1))/
434 : & (M(mu)*M(nu)*sqsin))-ReJq(mu,nu,qcount)
435 : ENDIF !invs
436 : ENDDO !nu
437 : mu=mu+nmagtype(imt)
438 : ENDDO !mu
439 :
440 : ENDIF !if(n.eq.1)
441 : ENDDO !n
442 : WRITE(oUnit,5006)enmax,qssmax
443 : WRITE(oUnit,5007)enmin,qssmin
444 : 5006 FORMAT('Maximal calculated energy = ',f20.10,'htr',
445 : & ' for the spin-spiral vector qss = ',3(f14.10,1x))
446 : 5007 FORMAT('Minimal calculated energy = ',f20.10,'htr',
447 : & ' for the spin-spiral vector qss = ',3(f14.10,1x))
448 :
449 : CLOSE(116)
450 : OPEN(117,file='shells',status='unknown')
451 : OPEN(118,file='MCinp',status='unknown')
452 : mu=1
453 :
454 : DO imt=1,mtypes
455 : DO nu=mu,nmagn
456 : WRITE(115,5004) itype(mu),itype(nu)
457 : WRITE(117,5004) itype(mu),itype(nu)
458 : 5004 FORMAT('#',i4,i4)
459 :
460 : sneq=0
461 : do ii=1,itype(mu)-1
462 : sneq=sneq+neq(ii)
463 : enddo
464 : if(itype(mu).le.sneq) then
465 : itype(mu)=sneq+1
466 : endif
467 :
468 : do ii=itype(mu),itype(nu)-1
469 : sneq=sneq+neq(ii)
470 : enddo
471 : if(itype(nu).le.sneq) then
472 : itype(nu)=sneq+1
473 : endif
474 :
475 : t(:)=taual(:,itype(mu))-taual(:,itype(nu))
476 :
477 : deltaz=taual(3,itype(mu))-taual(3,itype(nu)) !Added for film 07/10 S. Schroeder
478 : ! zcoord=taual(3,itype(nu))*amat(3,3)-taual(3,itype(1))*amat(3,3)
479 : zcoord=taual(3,itype(nu))*amat(3,3)
480 : tauC1(:)=taual(1,itype(mu))*amat(:,1)
481 : & +taual(2,itype(mu))*amat(:,2)
482 : & +taual(3,itype(mu))*amat(:,3)
483 : c... Generate the shells of neighbouring atoms
484 : CALL nshell(
485 : > amat,t,nsh,dims,nmax,shmax,film,zcoord,
486 : < nat,R,lenR,nop,mrot,deltaz)
487 :
488 : c ...
489 : c ... For the case of calculation with the least squares method
490 : c ... for one magnetic atom per unit cell
491 : IF (nshort.GT.0) THEN
492 : qcount=nqpt-1
493 : lwork=2*nshort
494 : ALLOCATE (Cmat(qcount,nshort),DelE(qcount),work(lwork))
495 : Cmat=0.0
496 : IF (nshort.GE.nqpt)THEN
497 : WRITE(*,*) ' Please supply the data for', nshort,
498 : & 'q-points different from zero'
499 : STOP
500 : ENDIF
501 :
502 : DO n=1,qcount
503 : DO nn=1,nshort
504 : DO atsh=1,nat(nn)
505 : scp=(q(1,1,n)*R(1,atsh,nn)
506 : & +q(2,1,n)*R(2,atsh,nn)
507 : & +q(3,1,n)*R(3,atsh,nn))*tpi
508 : Cmat(n,nn)=Cmat(n,nn)-1.0+cos(scp)
509 : ENDDO
510 : ENDDO
511 : DelE(n)=ReJq(1,1,n)*2000 ! multiply by 2000 to get [mRy/muB**2]
512 : ENDDO
513 :
514 : IF (nu.gt.mu) STOP 'For more than one magnetic atom in the unit
515 : & cell please set nshort=0'
516 : CALL dgels('N',qcount,nshort,1,Cmat,qcount,DelE,qcount,
517 : & work,lwork,info)
518 :
519 : c The routine dgels returns the solution, J(n), in the array DelE(n)
520 : Tc=0.0
521 : DO n=1,nshort
522 : Tc=Tc+nat(n)*DelE(n) !Mean-field Tc=1/3*(Sum_i(J_0,i))
523 : WRITE(115,5005) n,lenR(n),DelE(n) ! J in units [mRy/muB**2]
524 :
525 : c Now prepare and write the input for the Monte Carlo calculation:
526 : DO atsh=1,nat(n)
527 : tauC2(:)=(taual(1,itype(mu))-R(1,atsh,n))*amat(:,1)
528 : & +(taual(2,itype(mu))-R(2,atsh,n))*amat(:,2)
529 : & +(taual(3,itype(mu))-R(3,atsh,n))*amat(:,3)
530 :
531 : WRITE(118,5008) itype(mu),itype(nu),tauC1(1),tauC1(2),tauC1(3),
532 : & tauC2(1),tauC2(2),tauC2(3),DelE(n)
533 : ENDDO ! atsh
534 :
535 : ENDDO ! n
536 : DEALLOCATE (Cmat,DelE,work)
537 : ELSE
538 : c... Perform the back-Fourier transform
539 : DO nnn=1,nsh
540 : wrJ=0
541 : DO atsh=1,nat(nnn)
542 : IF(atsh.gt.shmax) STOP 'jcoff2:increase shmax!'
543 : J=0.0
544 : DO n=1,nqpt-1
545 : DO nn=1,nop
546 : IF(w(nn,n).EQ.1)THEN
547 : scp=(q(1,nn,n)*R(1,atsh,nnn)
548 : & +q(2,nn,n)*R(2,atsh,nnn)
549 : & +q(3,nn,n)*R(3,atsh,nnn))*tpi
550 : J=J+(((cos(scp))*ReJq(mu,nu,n)
551 : & +(sin(scp))*ImJq(mu,nu,n)))
552 : ENDIF
553 : ENDDO !nn
554 : ENDDO !n (qpts)
555 : J=(J/float(nqvect))*2000.0 ! J in units [mRy/muB**2]
556 : DO i=1,wrJ !A check for non-equivalent sub-shells
557 : IF(ABS(J-Jw(i)).LE.(tol))GOTO 55
558 : ENDDO
559 : WRITE(115,5005) nnn,lenR(nnn),J
560 : wrJ=wrJ+1
561 : Jw(wrJ)=J
562 : 55 CONTINUE
563 : 5005 FORMAT(i4,1x,2(f14.10,1x))
564 :
565 : c Now prepare and write the input for the Monte Carlo calculation:
566 : tauC2(:)=(taual(1,itype(mu))-R(1,atsh,nnn))*amat(:,1)
567 : & +(taual(2,itype(mu))-R(2,atsh,nnn))*amat(:,2)
568 : & +(taual(3,itype(mu))-R(3,atsh,nnn))*amat(:,3)
569 :
570 : WRITE(118,5008) itype(mu),itype(nu),tauC1(1),tauC1(2),tauC1(3),
571 : & tauC2(1),tauC2(2),tauC2(3),J
572 : ENDDO !atsh (atoms in each shell)
573 : c... In case of only one magnetic atom per unit cell, calculate the mean-field Tc
574 : IF(nmagn.EQ.1) Tc=Tc+nat(nnn)*J
575 : ENDDO !nnn (shells)
576 : ENDIF !nshort
577 : ENDDO !nu
578 : mu=mu+nmagtype(imt)
579 : ENDDO !imt
580 : IF(nmagn.EQ.1) THEN
581 : Tc=157.889*M(1)*M(1)*Tc/3.0
582 : WRITE(115,*) '# Tc(mean field)= ',Tc
583 : ENDIF
584 : 5008 FORMAT(i4,i4,7(1x,f14.10))
585 : CLOSE(117)
586 : CLOSE(118)
587 : CLOSE(115)
588 : #endif
589 0 : END SUBROUTINE priv_analyse_data
590 :
591 :
592 0 : END MODULE m_types_jij
|
