Commit a8abe8cf authored by Gregor Michalicek's avatar Gregor Michalicek

Added charge density I/O wrappers to some more routines

parent 5c45bce8
......@@ -63,7 +63,7 @@ MODULE m_cdn_io
! local variables
INTEGER :: mode, datend, k, i, iVac, j, iUnit
LOGICAL :: l_exist
LOGICAL :: l_exist, l_rhomatFile
CHARACTER(len=30) :: filename
CALL getMode(mode)
......@@ -96,9 +96,15 @@ MODULE m_cdn_io
IF (mode.EQ.CDN_DIRECT_MODE) THEN
filename = 'cdn1'
l_rhomatFile = .FALSE.
IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
INQUIRE(file="rhomat_inp",EXIST=l_exist)
IF (l_exist) filename = 'rhomat_inp'
IF (inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) THEN
INQUIRE(file="rhomat_out",EXIST=l_exist)
IF (l_exist) filename = 'rhomat_out'
END IF
IF(l_exist) l_rhomatFile = .TRUE.
END IF
IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
filename = 'cdn'
......@@ -112,10 +118,7 @@ MODULE m_cdn_io
iUnit = 93
OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='old')
IF (inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) THEN
IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
CALL juDFT_error("inOrOutCDN.EQ.CDN_OUTPUT_DEN_const incompatible to l_noco",calledby ="readDensity")
END IF
IF ((inOrOutCDN.EQ.CDN_OUTPUT_DEN_const).AND.(archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const)) THEN
! call loddop to move the file position to the output density
CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
iUnit,iter,fr,fpw,fz,fzxy)
......@@ -126,7 +129,7 @@ MODULE m_cdn_io
iUnit,iter,fr,fpw,fz,fzxy)
! read in additional data if l_noco and data is present
IF ((archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const).AND.(TRIM(ADJUSTL(filename)).EQ.'rhomat_inp')) THEN
IF ((archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const).AND.l_rhomatFile) THEN
READ (iUnit,iostat=datend) (cdom(k),k=1,stars%ng3)
IF (datend == 0) THEN
IF (input%film) THEN
......@@ -204,6 +207,7 @@ MODULE m_cdn_io
filename = 'cdn1'
IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
filename = 'rhomat_inp'
IF(inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) filename = 'rhomat_out'
END IF
IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
filename = 'cdn'
......@@ -212,10 +216,7 @@ MODULE m_cdn_io
iUnit = 93
OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='unknown')
IF (inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) THEN
IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
CALL juDFT_error("inOrOutCDN.EQ.CDN_OUTPUT_DEN_const incompatible to l_noco",calledby ="writeDensity")
END IF
IF ((inOrOutCDN.EQ.CDN_OUTPUT_DEN_const).AND.(archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const)) THEN
! Generate data in temp arrays and variables to be able to perform loddop call.
! loddop is called to move the file position to the output density position.
......
......@@ -21,6 +21,7 @@
USE m_prpqfftmap
USE m_cdnval
USE m_loddop
USE m_cdn_io
USE m_wrtdop
USE m_cdntot
USE m_cdnovlp
......@@ -59,7 +60,7 @@
REAL fix,qtot,scor,seig,smom,stot,sval,dummy
REAL slmom,slxmom,slymom,sum,thetai,phii
INTEGER iter,ivac,j,jspin,jspmax,k,n,nt,ieig,ikpt
INTEGER ityp,ilayer,urec,itype,iatom
INTEGER ityp,ilayer,urec,itype,iatom,archiveType
LOGICAL l_relax_any,exst,n_exist,l_st
TYPE(t_noco)::noco_new
! ..
......@@ -101,13 +102,12 @@
!
! Read in input density
!
nt = 71
IF ( (.NOT. noco%l_noco) .AND. mpi%irank.EQ.0) THEN
OPEN (nt,file='cdn1',form='unformatted',status='old')
REWIND (nt)
CALL loddop(stars,vacuum,atoms,sphhar, input,sym, nt, iter,rho,qpw,rht,rhtxy)
ENDIF
!
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
IF(noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
ALLOCATE(cdom(1),cdomvz(1,1),cdomvxy(1,1,1))
CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
CDN_INPUT_DEN_const,0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
DEALLOCATE(cdom,cdomvz,cdomvxy)
IF (mpi%irank.EQ.0) THEN
INQUIRE(file='enpara',exist=l_enpara)
......@@ -473,27 +473,9 @@ enddo
CLOSE(20)
CALL juDFT_error("slice OK",calledby="cdngen")
END IF
!
IF (noco%l_noco) THEN
!---> pk non-collinear
!---> write output density matrix on file rhomat_out
!---> first the diagonal elements of the density matrix
OPEN (26,FILE='rhomat_out',FORM='unformatted',&
& STATUS='unknown')
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 26, iter,rho,qpw,rht,rhtxy)
!---> and then the off-diagonal part
WRITE (26) (cdom(k),k=1,stars%ng3)
IF (input%film) THEN
WRITE (26) ((cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac)
WRITE (26) (((cdomvxy(j,k-1,ivac),j=1,vacuum%nmzxy), k=2,oneD%odi%nq2),ivac=1,vacuum%nvac)
ENDIF
CLOSE (26)
!---> pk non-collinear
ELSE
! ----> write output density on unit 71
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, nt, iter,rho,qpw,rht,rhtxy)
CLOSE (nt)
ENDIF
CALL writeDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
CDN_OUTPUT_DEN_const,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
ENDIF
DEALLOCATE (cdom,cdomvz,cdomvxy,qa21)
......
......@@ -16,7 +16,7 @@ CONTAINS
SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym, cell, it, noco, oneD,hybrid)
!
#include"cpp_double.h"
USE m_loddop
USE m_cdn_io
USE m_wrtdop
USE m_brysh1
USE m_stmix
......@@ -49,7 +49,7 @@ CONTAINS
REAL fix,intfac,vacfac
INTEGER i,iter,imap,js,mit,nt,nt1,irecl
INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
INTEGER iq2,iq3,ivac,imz ,iofl
INTEGER iq2,iq3,ivac,imz ,iofl, archiveType
INTEGER n_u_keep
LOGICAL lexist,l_ldaU
INTEGER d1,d10,asciioffset
......@@ -143,9 +143,18 @@ CONTAINS
IF (noco%l_noco) THEN
ALLOCATE (cdom(stars%n3d),cdomvz(vacuum%nmzd,2), cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2))
archiveType = CDN_ARCHIVE_TYPE_NOCO_const
ELSE
ALLOCATE ( cdom(1),cdomvz(1,1),cdomvxy(1,1,1) )
ALLOCATE (cdom(1),cdomvz(1,1),cdomvxy(1,1,1))
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
ENDIF
!---> initialize arrays for the off-diagonal part of the density matrix
cdom(:) = CMPLX(0.0,0.0)
IF (input%film) THEN
cdomvz(:,:) = CMPLX(0.0,0.0)
cdomvxy(:,:,:) = CMPLX(0.0,0.0)
END IF
!
INQUIRE (file='broyd.'//CHAR(input%imix+48),exist=lexist)
DO i = 1,6
......@@ -193,44 +202,15 @@ CONTAINS
!gs-
!---> reload densities of current iteration
IF (noco%l_noco) THEN
!---> initialize arrays for the off-diagonal part of the density matrix
cdom(:) = CMPLX(0.0,0.0)
IF (input%film) THEN
cdomvz(:,:) = CMPLX(0.0,0.0)
cdomvxy(:,:,:) = CMPLX(0.0,0.0)
END IF
CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
CDN_INPUT_DEN_const,0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
!---> reload input density matrix from file rhomat_inp
OPEN (nrhomfile,FILE='rhomat_inp',FORM='unformatted', STATUS='unknown')
!---> first the diagonal elements of the density matrix
CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
nrhomfile, iter,rho,qpw,rht,rhtxy)
!---> and then the off-diagonal part
READ (nrhomfile,END=150,ERR=50) (cdom(iq3),iq3=1,stars%ng3)
IF (input%film) THEN
READ (nrhomfile,END=75,ERR=50) ((cdomvz(imz,ivac), imz=1,vacuum%nmz),ivac=1,vacuum%nvac)
READ (nrhomfile,END=75,ERR=50) (((cdomvxy(imz,iq2-1,ivac), imz=1,vacuum%nmzxy),iq2=2,oneD%odi%nq2),ivac=1,vacuum%nvac)
ENDIF
GOTO 150
50 WRITE(6,*)'rhodirgen: ERROR: Problems while reading density'
WRITE(6,*)'matrix from file rhomat_inp.'
CALL juDFT_error("ERROR while reading file rhomat_inp",calledby ="mix")
75 WRITE(6,*)'rhomatdir: ERROR: reached end of file rhomat_inp'
WRITE(6,*)'while reading the vacuum part of the off-diagonal'
WRITE(6,*)'element of the desity matrix.'
CALL juDFT_error("ERROR while reading file rhomat_inp",calledby ="mix")
150 CLOSE (nrhomfile)
ELSE
nt = 71 !gs see above
OPEN (nt,file='cdn1',form='unformatted',status='old')
REWIND nt
CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
nt, iter,rho,qpw,rht,rhtxy)
IF (.NOT.noco%l_noco) THEN
!---> write density to file for storage
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
nt1, iter,rho,qpw,rht,rhtxy)
ENDIF
END IF
!
!---> put input charge density into arrays sm
! in the spin polarized case the arrays consist of
......@@ -238,29 +218,17 @@ CONTAINS
CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
intfac,vacfac,qpw,rho,rht,rhtxy,cdom,cdomvz,cdomvxy,n_mmp(-3,-3,1,1,1), nmap,nmaph,mapmt,mapvac,mapvac2,sm)
! load output charge density
!
IF (noco%l_noco) THEN
!---> reload output density matrix from file rhomat_out
OPEN (nrhomfile,FILE='rhomat_out',FORM='unformatted', STATUS='unknown')
!---> first the diagonal elements of the density matrix
CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
nrhomfile, iter,rho,qpw,rht,rhtxy)
!---> and then the off-diagonal part
READ (nrhomfile) (cdom(iq3),iq3=1,stars%ng3)
IF (input%film) THEN
READ (nrhomfile) ((cdomvz(imz,ivac),imz=1,vacuum%nmz),ivac=1,vacuum%nvac)
READ (nrhomfile) (((cdomvxy(imz,iq2-1,ivac),imz=1,vacuum%nmzxy), iq2=2,oneD%odi%nq2),ivac=1,vacuum%nvac)
ENDIF
CLOSE (nrhomfile)
ELSE
CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
nt, iter,rho,qpw,rht,rhtxy)
CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
CDN_OUTPUT_DEN_const,0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
IF (.NOT.noco%l_noco) THEN
!---> write density to file for storage
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
nt1, iter,rho,qpw,rht,rhtxy)
CLOSE(nt1)
ENDIF
END IF
!
!---> put output charge density into arrays fsm
!
......@@ -367,31 +335,10 @@ CONTAINS
qpw,rhtxy,rho,rht,.FALSE., fix)
iter = iter + 1
IF (noco%l_noco) THEN
!
!---> write mixed density to file rhomat_out
!
OPEN (nrhomfile,FILE='rhomat_inp',FORM='unformatted', STATUS='unknown')
!---> first the diagonal elements of the density matrix
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
nrhomfile, iter,rho,qpw,rht,rhtxy)
!---> and then the off-diagonal part
WRITE (nrhomfile) (cdom(iq3),iq3=1,stars%ng3)
IF (input%film) THEN
WRITE (nrhomfile) ((cdomvz(imz,ivac),imz=1,vacuum%nmz),ivac=1,vacuum%nvac)
WRITE (nrhomfile) (((cdomvxy(imz,iq2-1,ivac),imz=1,vacuum%nmzxy),&
iq2=2,oneD%odi%nq2),ivac=1,vacuum%nvac)
ENDIF
CLOSE (nrhomfile)
ELSE
!
!---> write new density onto unit 71 (overwrite)
!
REWIND nt
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
nt, iter,rho,qpw,rht,rhtxy)
CLOSE (nt)
ENDIF
CALL writeDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
DEALLOCATE ( cdom,cdomvz,cdomvxy )
IF ( atoms%n_u > 0 ) THEN
OPEN (69,file='n_mmp_mat',status='replace',form='formatted')
......
......@@ -169,9 +169,7 @@ CONTAINS
CALL timestart("optional: spin polarized density")
CALL cdnsp(&
& atoms,input,vacuum,sphhar,&
& stars,&
& sym,&
& cell,DIMENSION)
& stars,sym,oneD,cell,dimension)
!
CALL timestop("optional: spin polarized density")
END IF
......@@ -211,7 +209,7 @@ CONTAINS
IF (input%l_bmt) THEN
CALL bmt(&
& stars,input,noco,atoms,sphhar,vacuum,&
& cell,sym)
& cell,sym,oneD)
ENDIF
ENDIF ! mpi%irank == 0
......
......@@ -46,6 +46,7 @@ CONTAINS
USE m_force_a3
USE m_forcew
USE m_loddop
USE m_cdn_io
USE m_icorrkeys
USE m_types
USE m_xmlOutput
......@@ -69,11 +70,12 @@ CONTAINS
! ..
! .. Local Scalars ..
REAL rhs,totz, eigSum
INTEGER n,j,nt,iter,i
INTEGER n,j,nt,iter,i, archiveType
! .. Local Arrays ..
REAL vmd(atoms%ntype),zintn_r(atoms%ntype)
REAL dpj(atoms%jmtd)
COMPLEX :: cdom(1),cdomvz(1,1),cdomvxy(1,1,1)
CHARACTER(LEN=20) :: attributes(3)
!.....density
REAL, ALLOCATABLE :: rho(:,:,:,:),rht(:,:,:)
......@@ -139,16 +141,12 @@ CONTAINS
! ----> VM terms
! ---> reload the density
!
IF (noco%l_noco) THEN
nt = 70
OPEN (nt,file='cdn',form='unformatted',status='old')
ELSE
nt = 71
OPEN (nt,file='cdn1',form='unformatted',status='old')
ENDIF
CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
nt, iter,rho,qpw,rht,rhtxy)
CLOSE (nt)
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_CDN_const
CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
CDN_INPUT_DEN_const,0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
!+for
! ---> reload the COULOMB potential
!
......
......@@ -2,11 +2,11 @@ MODULE m_bmt
contains
SUBROUTINE bmt(&
& stars,input,noco,atoms,sphhar,vacuum,&
& cell,sym)
& cell,sym,oneD)
!
use m_types
use m_juDFT
USE m_loddop
USE m_cdn_io
USE m_wrtdop
IMPLICIT NONE
! ..
......@@ -18,11 +18,14 @@ contains
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_oneD),INTENT(IN) :: oneD
INTEGER k,i,ivac ,it
INTEGER type,typmag
INTEGER type,typmag, archiveType
CHARACTER(len=8) filename
COMPLEX, ALLOCATABLE :: fpw(:,:),fzxy(:,:,:,:)
REAL, ALLOCATABLE :: fz(:,:,:),fr(:,:,:,:)
COMPLEX, ALLOCATABLE :: cdom(:),cdomvz(:,:),cdomvxy(:,:,:)
! ..
! ..
......@@ -38,19 +41,15 @@ contains
!atoms%jmtd = maxval(atoms%jri(:))
!sphhar%nlhd = maxval(sphhar%nlh(:))
ALLOCATE( fpw(stars%ng3,input%jspins),fzxy(vacuum%nmzxy,stars%ng2-1,2,input%jspins) )
ALLOCATE( fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins),fz(vacuum%nmz,2,input%jspins) )
ALLOCATE(fpw(stars%ng3,input%jspins),fzxy(vacuum%nmzxy,stars%ng2-1,2,input%jspins))
ALLOCATE(fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins),fz(vacuum%nmz,2,input%jspins))
ALLOCATE(cdom(stars%n3d),cdomvz(vacuum%nmzd,2),cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2))
IF (noco%l_noco) THEN
OPEN (98,file='rhomat_inp',form='unformatted',action='read')
ELSE
OPEN (98,file='cdn1',form='unformatted',action='read')
ENDIF
CALL loddop(&
& stars,vacuum,atoms,sphhar,input,sym,&
& 98,&
& it,fr,fpw,fz,fzxy)
CLOSE (98)
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
CDN_INPUT_DEN_const,0,it,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
IF ( typmag < atoms%ntype ) THEN
DO type= typmag+1,atoms%ntype
......@@ -97,7 +96,8 @@ contains
& it,fr,fpw,fz,fzxy)
CLOSE(98)
DEALLOCATE( fpw,fzxy,fr,fz)
DEALLOCATE(cdom,cdomvz,cdomvxy)
DEALLOCATE(fpw,fzxy,fr,fz)
END SUBROUTINE bmt
END MODULE m_bmt
......@@ -15,12 +15,11 @@
CONTAINS
SUBROUTINE cdnsp(&
& atoms,input,vacuum,sphhar,&
& stars,sym,cell,DIMENSION)
& stars,sym,oneD,cell,DIMENSION)
USE m_intgr, ONLY : intgr3
USE m_constants, ONLY : pi_const
USE m_loddop
USE m_wrtdop
USE m_cdn_io
USE m_types
IMPLICIT NONE
! ..
......@@ -30,17 +29,19 @@
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_input),INTENT(INOUT) :: input
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_dimension),INTENT(IN)::DIMENSION
TYPE(t_dimension),INTENT(IN) :: DIMENSION
! ..
! .. Local Scalars ..
REAL dummy,p,pp,qtot1,qtot2,spmtot,qval,sfp
INTEGER i,iter,ivac,j,k,lh,n,na,nt,jsp_new
INTEGER i,iter,ivac,j,k,lh,n,na,jsp_new
INTEGER ios
LOGICAL n_exist
! ..
! .. Local Arrays ..
REAL rhoc(DIMENSION%msh,atoms%ntype)
COMPLEX :: cdom(1),cdomvz(1,1),cdomvxy(1,1,1)
COMPLEX, ALLOCATABLE :: qpw(:,:),rhtxy(:,:,:,:)
REAL , ALLOCATABLE :: rho(:,:,:,:),rht(:,:,:)
CHARACTER(len=140), ALLOCATABLE :: clines(:)
......@@ -63,19 +64,11 @@
ENDDO
CLOSE (17)
nt = 71 !gs, see sub mix
OPEN (nt,file='cdn1',form='unformatted',status='old')
!
! ---> set jspins=1 to read the paramagnetic density
!
input%jspins=1
CALL loddop(&
& stars,vacuum,atoms,sphhar,&
& input,sym,&
& nt,&
& iter,rho,qpw,rht,rhtxy)
CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
CDN_INPUT_DEN_const,0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
input%jspins=2
!
qval = 0.
na = 1
!
......@@ -120,14 +113,9 @@
ENDDO
ENDDO
ENDIF
! ----> write the spin-polarized density on unit nt
REWIND nt
CALL wrtdop(&
& stars,vacuum,atoms,sphhar,&
& input,sym,&
& nt,&
& iter,rho,qpw,rht,rhtxy)
CLOSE (nt)
! ----> write the spin-polarized density
CALL writeDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
CDN_INPUT_DEN_const,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
!
! -----> This part is only used for testing th e magnetic moment in
! -----> each sphere
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment