Commit 248f92eb authored by Daniel Wortmann's avatar Daniel Wortmann

Added code for fixing the density after shifting the atoms

parent afbbfd30
......@@ -17,6 +17,7 @@ force/force_a8.F90
force/force_b8.f90
force/force_sf.F90
force/force_w.F90
force/fix_by_gaussian.F90
#force/geo.f90
force/stern.f90
force/relaxation.F90
......
!--------------------------------------------------------------------------------
! Copyright (c) 2019 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_fix_by_gaussian
USE m_judft
IMPLICIT NONE
CONTAINS
SUBROUTINE fix_by_gaussian(shift,atoms,stars,mpi,sym,vacuum,sphhar,input,oned,cell,noco,den)
!The idea of this fix is to add an Gaussian to the INT which make the charge flat at the MT-boundary and to
!shift this Gaussian with the displacement.
USE m_qfix
USE m_spgrot
USE m_constants
USE m_types
IMPLICIT NONE
REAL,INTENT(in) :: shift(:,:)
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_input),INTENT(IN) :: input
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_potden),INTENT(INOUT):: den
REAL :: slope1,slope2,dr,alpha
REAL :: sigma,a_fac,gauss,x,fix
INTEGER :: kr(3,sym%nop)
COMPLEX :: sf,phas(sym%nop)
INTEGER :: js,n,l,k,nat,j
TYPE(t_input) :: inputtmp
DO js=1,input%jspins
DO n=1,atoms%ntype
DO l=0,0 !Currently only l=0
! alpha = LOG( den%mt(atoms%jri(n)-1,l,n,js) / den%mt(atoms%jri(n),l,n,js) )
! alpha = SQRT(alpha / ( atoms%rmt(n)*atoms%rmt(n)*( 1.0-EXP( -2.0*atoms%dx(n) ) ) ))
! A_fac= den%mt(atoms%jri(n),l,n,js)/gaussian_r(atoms%rmt(n),alpha)
dr=atoms%rmsh(atoms%jri(n)-1,n)-atoms%rmsh(atoms%jri(n),n)
slope1=(den%mt(atoms%jri(n)-1,l,n,js)/atoms%rmsh(atoms%jri(n)-1,n)**2-den%mt(atoms%jri(n),l,n,js)/atoms%rmsh(atoms%jri(n),n)**2)/dr
slope1=(den%mt(atoms%jri(n)-1,l,n,js)-den%mt(atoms%jri(n),l,n,js))/dr/atoms%rmt(n)**2
sigma=atoms%rmt(n)/2.0
alpha=1/sigma
slope2=(gaussian_r(atoms%rmsh(atoms%jri(n)-1,n),alpha)-gaussian_r(atoms%rmsh(atoms%jri(n),n),alpha))/dr
A_fac=slope1/slope2
PRINT *, a_fac, 1/alpha
DO k=2,stars%ng3
gauss=A_fac*gaussian_g(stars%sk3(k),alpha)
CALL spgrot(sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,stars%kv3(:,k),kr,phas)
DO nat=SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
sf=0.0
DO j = 1,sym%nop
x=-tpi_const*DOT_PRODUCT(1.*kr(:,j),atoms%taual(:,nat))
sf = sf + CMPLX(COS(x),SIN(x))*CONJG(phas(j))
x=-tpi_const*DOT_PRODUCT(1.*kr(:,j),atoms%taual(:,nat)+shift(:,nat))
sf = sf - CMPLX(COS(x),SIN(x))*CONJG(phas(j))
ENDDO
ENDDO
den%pw(k,js)=den%pw(k,js)+gauss*sf/sym%nop/cell%omtil
ENDDO
ENDDO
END DO
END DO
inputtmp=input
inputtmp%qfix=1
CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,inputtmp,cell,oneD,&
den,noco%l_noco,mpi%isize==1,force_fix=.TRUE.,fix=fix)
END SUBROUTINE fix_by_gaussian
FUNCTION gaussian_r(r,alpha)
REAL,INTENT(IN) :: r,alpha
real :: gaussian_r
gaussian_r=EXP(-r**2*alpha**2)
END FUNCTION gaussian_r
FUNCTION gaussian_g(g,alpha)
USE m_constants
REAL,INTENT(IN) :: g,alpha
REAL :: gaussian_g
gaussian_g=SQRT(pi_const**3)/alpha**3*EXP(-0.25*g**2/alpha**2)
END FUNCTION gaussian_g
END MODULE m_fix_by_gaussian
......@@ -44,10 +44,10 @@ CONTAINS
if (l_firstcall) THEN
INQUIRE(file='qfix',exist=l_qfixfile)
IF (l_qfixfile) CALL judft_info("qfix file no longer supported, check the qfix option in inp.xml","INFO")
if (input%qfix==1) fixtotal=.false.
IF (input%qfix==1) fixtotal=.FALSE.
l_firstcall=.false.
ELSE
IF (input%qfix==0.AND..NOT.force_fix) RETURN
IF (MOD(input%qfix,2)==0.AND..NOT.force_fix) RETURN
ENDIF
! qfix==0 means no qfix was given in inp.xml.
! In this case do nothing except when forced to fix!
......
......@@ -142,7 +142,7 @@
! Store structure data
CALL storeStructureIfNew(input, atoms, cell, vacuum, oneD, sym)
CALL storeStructureIfNew(input,stars, atoms, cell, vacuum, oneD, sym,mpi,sphhar,noco)
!+odim
IF (input%film.OR.(sym%namgrp.NE.'any ')) THEN
......
......@@ -502,7 +502,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
! Store structure data
CALL storeStructureIfNew(input, atoms, cell, vacuum, oneD, sym)
CALL storeStructureIfNew(input,stars, atoms, cell, vacuum, oneD, sym, mpi,sphhar,noco)
! Generate stars
......@@ -553,6 +553,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
END IF
CALL transform_by_moving_atoms(mpi,stars,atoms,vacuum, cell, sym, sphhar,input,oned,noco)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! End of input postprocessing (calculate missing parameters)
......
......@@ -34,7 +34,7 @@ MODULE m_cdn_io
PUBLIC readStars, writeStars
PUBLIC readStepfunction, writeStepfunction
PUBLIC setStartingDensity, readPrevEFermi, deleteDensities
PUBLIC storeStructureIfNew
PUBLIC storeStructureIfNew,transform_by_moving_atoms
PUBLIC getIOMode
PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
......@@ -841,7 +841,7 @@ MODULE m_cdn_io
END SUBROUTINE writeCoreDensity
SUBROUTINE storeStructureIfNew(input, atoms, cell, vacuum, oneD, sym)
SUBROUTINE storeStructureIfNew(input,stars, atoms, cell, vacuum, oneD, sym,mpi,sphhar,noco)
TYPE(t_input),INTENT(IN) :: input
TYPE(t_atoms), INTENT(IN) :: atoms
......@@ -849,6 +849,10 @@ MODULE m_cdn_io
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_input) :: inputTemp
TYPE(t_atoms) :: atomsTemp
......@@ -885,6 +889,7 @@ MODULE m_cdn_io
END IF
END IF
IF (l_writeStructure) THEN
CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, sym, currentStructureIndex,.TRUE.)
CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
......@@ -902,6 +907,105 @@ MODULE m_cdn_io
END SUBROUTINE storeStructureIfNew
SUBROUTINE transform_by_moving_atoms(mpi,stars,atoms,vacuum, cell, sym, sphhar,input,oned,noco)
USE m_types
USE m_qfix
use m_fix_by_gaussian
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_input),INTENT(IN) :: input
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
!Locals
INTEGER :: archiveType
LOGICAL :: l_qfix
REAL :: fermiEnergy,fix
REAL :: shifts(3,SIZE(atoms%taual,2))
TYPE(t_potden):: den
TYPE(t_input) :: inputTemp
TYPE(t_atoms) :: atomsTemp
TYPE(t_cell) :: cellTemp
TYPE(t_vacuum) :: vacuumTemp
TYPE(t_oneD) :: oneDTemp
TYPE(t_sym) :: symTemp
INTEGER :: mode,currentStarsIndex,currentLatharmsIndex,currentStructureIndex
INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex,structureindex
LOGICAL :: l_same,l_structure_by_shift
#ifdef CPP_HDF
INTEGER(HID_T) :: fileID
character(len=50) :: archivename
#endif
#ifdef CPP_MPI
INCLUDE 'mpif.h'
integer :: ierr
#endif
l_same=.TRUE.;l_structure_by_shift=.TRUE.
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
IF (mpi%irank==0) THEN
CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
IF (currentStructureIndex>0.AND.lastdensityindex>0) THEN
!Determine structure of last density
WRITE(archivename,"(a,i0)") "cdn-",lastdensityindex
CALL peekDensityEntryHDF(fileID, archivename, DENSITY_TYPE_IN_const, structureIndex=structureIndex)
!Read that structure
CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, symTemp, StructureIndex)
CALL compareStructure(input, atoms, vacuum, cell, sym, inputTemp, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same,l_structure_by_shift)
ENDIF
CALL closeCDNPOT_HDF(fileID)
ENDIF
#ifdef CPP_MPI
CALL mpi_bc(l_same,1,MPI_LOGICAL,0,ierr)
CALL mpi_bc(l_structure_by_shift,1,MPI_LOGICAL,0,ierr)
#endif
IF (l_same.OR..NOT.l_structure_by_shift) RETURN ! nothing to do
IF (mpi%irank==0) THEN
WRITE(6,*) "Atomic movement detected, trying to adjust charge density"
!Calculate shifts
shifts=atomsTemp%taual-atoms%taual
!Determine type of charge
archiveType = MERGE(CDN_ARCHIVE_TYPE_NOCO_const,CDN_ARCHIVE_TYPE_CDN1_const,noco%l_noco)
!read the current density
CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,fermiEnergy,l_qfix,den)
ENDIF
!Now fix the density
SELECT CASE(input%qfix)
CASE (0,1) !just qfix the density
IF (mpi%irank==0) WRITE(6,*) "Using qfix to adjust density"
CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
den,noco%l_noco,mpi%isize==1,force_fix=.TRUE.,fix=fix)
CASE(2)
CALL fix_by_gaussian(shifts,atoms,stars,mpi,sym,vacuum,sphhar,input,oned,cell,noco,den)
CASE default
CALL judft_error("Wrong choice of qfix in input")
END SELECT
!Now write the density to file
IF (mpi%irank==0) CALL writedensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,-1.0,fermiEnergy,l_qfix,den)
#endif
END IF
END SUBROUTINE transform_by_moving_atoms
SUBROUTINE writeStars(stars,l_xcExtended,l_ExtData)
TYPE(t_stars),INTENT(IN) :: stars
......
......@@ -61,7 +61,7 @@ MODULE m_cdnpot_io_common
END SUBROUTINE compareStepfunctions
SUBROUTINE compareStructure(input, atoms, vacuum, cell, sym, refInput, refAtoms, refVacuum,&
refCell, refSym, l_same)
refCell, refSym, l_same,l_shift_only)
TYPE(t_input),INTENT(IN) :: input, refInput
TYPE(t_atoms),INTENT(IN) :: atoms, refAtoms
......@@ -70,10 +70,12 @@ MODULE m_cdnpot_io_common
TYPE(t_sym),INTENT(IN) :: sym, refSym
LOGICAL, INTENT(OUT) :: l_same
LOGICAL,OPTIONAL,INTENT(OUT) ::l_shift_only
INTEGER :: i
l_same = .TRUE.
IF(atoms%ntype.NE.refAtoms%ntype) l_same = .FALSE.
IF(atoms%nat.NE.refAtoms%nat) l_same = .FALSE.
......@@ -97,6 +99,9 @@ MODULE m_cdnpot_io_common
IF(ANY(sym%mrot(:,:,:sym%nop).NE.refSym%mrot(:,:,:sym%nop))) l_same = .FALSE.
IF(ANY(ABS(sym%tau(:,:sym%nop)-refSym%tau(:,:sym%nop)).GT.1e-10)) l_same = .FALSE.
END IF
IF (PRESENT(l_shift_only)) l_shift_only=l_same
!Now the positions are checked...
IF(l_same) THEN
DO i = 1, atoms%nat
IF(ANY(ABS(atoms%pos(:,i)-refAtoms%pos(:,i)).GT.1e-10)) l_same = .FALSE.
......
......@@ -2650,11 +2650,11 @@ MODULE m_cdnpot_io_hdf
INTEGER, INTENT(IN) :: densityType
CHARACTER(LEN=*), INTENT(IN) :: archiveName
INTEGER, INTENT(OUT) :: date, time, iter
INTEGER, INTENT(OUT) :: starsIndex, latharmsIndex, structureIndex, stepfunctionIndex
INTEGER, INTENT(OUT) :: previousDensityIndex, jspins
REAL, INTENT(OUT) :: fermiEnergy, distance
LOGICAL, INTENT(OUT) :: l_qfix
INTEGER, INTENT(OUT),OPTIONAL :: date, time, iter
INTEGER, INTENT(OUT),OPTIONAL :: starsIndex, latharmsIndex, structureIndex, stepfunctionIndex
INTEGER, INTENT(OUT),OPTIONAL :: previousDensityIndex, jspins
REAL, INTENT(OUT),OPTIONAL :: fermiEnergy, distance
LOGICAL, INTENT(OUT),OPTIONAL :: l_qfix
INTEGER :: localDensityType
LOGICAL :: l_exist
......@@ -2706,20 +2706,20 @@ MODULE m_cdnpot_io_hdf
CALL h5gopen_f(fileID, TRIM(ADJUSTL(archiveName)), archiveID, hdfError)
CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError)
CALL io_read_attint0(archiveID,'previousDensityIndex',previousDensityIndex)
CALL io_read_attint0(archiveID,'starsIndex',starsIndex)
CALL io_read_attint0(archiveID,'latharmsIndex',latharmsIndex)
CALL io_read_attint0(archiveID,'structureIndex',structureIndex)
CALL io_read_attint0(archiveID,'stepfunctionIndex',stepfunctionIndex)
CALL io_read_attint0(archiveID,'spins',jspins)
CALL io_read_attint0(archiveID,'iter',iter)
CALL io_read_attint0(archiveID,'date',date)
CALL io_read_attint0(archiveID,'time',time)
CALL io_read_attreal0(archiveID,'distance',distance)
IF (PRESENT(previousDensityIndex)) CALL io_read_attint0(archiveID,'previousDensityIndex',previousDensityIndex)
IF (PRESENT(starsIndex)) CALL io_read_attint0(archiveID,'starsIndex',starsIndex)
IF (PRESENT(latharmsIndex)) CALL io_read_attint0(archiveID,'latharmsIndex',latharmsIndex)
IF (PRESENT(structureIndex)) CALL io_read_attint0(archiveID,'structureIndex',structureIndex)
IF (PRESENT(stepfunctionIndex)) CALL io_read_attint0(archiveID,'stepfunctionIndex',stepfunctionIndex)
IF (PRESENT(jspins)) CALL io_read_attint0(archiveID,'spins',jspins)
IF (PRESENT(iter)) CALL io_read_attint0(archiveID,'iter',iter)
IF (PRESENT(date)) CALL io_read_attint0(archiveID,'date',date)
IF (PRESENT(time)) CALL io_read_attint0(archiveID,'time',time)
IF (PRESENT(distance)) CALL io_read_attreal0(archiveID,'distance',distance)
IF (densityType.NE.DENSITY_TYPE_UNDEFINED_const) THEN
CALL io_read_attreal0(groupID,'fermiEnergy',fermiEnergy)
CALL io_read_attlog0(groupID,'l_qfix',l_qfix)
IF (PRESENT(fermiEnergy)) CALL io_read_attreal0(groupID,'fermiEnergy',fermiEnergy)
IF (PRESENT(l_qfix)) CALL io_read_attlog0(groupID,'l_qfix',l_qfix)
END IF
CALL h5gclose_f(groupID, hdfError)
......
......@@ -671,14 +671,7 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
input%epsforce = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsforce'))
input%forcemix = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@forcemix'))
input%force_converged = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@force_converged'))
numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@qfix')
IF (numberNodes.EQ.1) THEN
input%qfix = 1
l_qfix = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@qfix'))
IF (l_qfix) THEN
input%qfix = 2
END IF
END IF
input%qfix = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@qfix'))
END IF
! Read in optional general LDA+U parameters
......
......@@ -590,8 +590,9 @@
<xsd:attribute name="forcealpha" type="xsd:string" use="required"/>
<xsd:attribute name="epsdisp" type="xsd:string" use="optional" default="0.001"/>
<xsd:attribute name="epsforce" type="xsd:string" use="optional" default="0.001"/>
<xsd:attribute name="forcemix" type="xsd:string" use="optional" default="3" />
<xsd:attribute name="qfix" type="FleurBool" use="optional"/>
<xsd:attribute name="forcemix" type="xsd:integer" use="optional" default="3" />
<xsd:attribute name="qfix" type="xsd:integer" use="optional" default="0" />
<xsd:attribute name="force_converged" type="xsd:string" use="optional" default="0.00001" />
</xsd:complexType>
<xsd:complexType name="LdaUGeneralType">
......
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