LCOV - code coverage report
Current view: top level - forcetheorem - jij.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 75 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 9 0.0 %

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

Generated by: LCOV version 1.13