Commit 13adb7a0 authored by Miriam Hinzen's avatar Miriam Hinzen

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents 2aa60e57 5175e5ca
......@@ -29,7 +29,7 @@ CONTAINS
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_potden),INTENT(IN) :: v
CLASS(t_mat),INTENT(INOUT) :: hmat(:,:),smat(:,:)
CLASS(t_mat),INTENT(INOUT) :: hmat(:,:),smat(:,:)
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp
......@@ -40,14 +40,14 @@ CONTAINS
! .. Local Scalars ..
COMPLEX hij,sij,apw_lo,c_1
REAL d2,gz,sign,th,wronk
INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ipot,ii0,i0
INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ii0,i0
INTEGER ivac,irec,imz,igvm2,igvm2i,s1,s2
INTEGER jspin1,jspin2,jmax,jsp_start,jsp_end
INTEGER i_start,nc,nc_0
! ..
! .. Local Arrays ..
INTEGER:: nv2(input%jspins)
INTEGER kvac1(DIMENSION%nv2d,input%jspins),kvac2(DIMENSION%nv2d,input%jspins)
INTEGER kvac(2,DIMENSION%nv2d,input%jspins)
INTEGER map2(DIMENSION%nvd,input%jspins)
COMPLEX tddv(DIMENSION%nv2d,DIMENSION%nv2d),tduv(DIMENSION%nv2d,DIMENSION%nv2d)
COMPLEX tudv(DIMENSION%nv2d,DIMENSION%nv2d),tuuv(DIMENSION%nv2d,DIMENSION%nv2d)
......@@ -66,16 +66,14 @@ CONTAINS
nv2(jspin) = 0
k_loop:DO k = 1,lapw%nv(jspin)
DO j = 1,nv2(jspin)
IF (lapw%k1(k,jspin).EQ.kvac1(j,jspin)&
.AND. lapw%k2(k,jspin).EQ.kvac2(j,jspin)) THEN
IF (all(lapw%gvec(1:2,k,jspin)==kvac(1:2,j,jspin))) THEN
map2(k,jspin) = j
CYCLE k_loop
END IF
ENDDO
nv2(jspin) = nv2(jspin) + 1
IF (nv2(jspin)>DIMENSION%nv2d) CALL juDFT_error("hsvac:dimension%nv2d",calledby ="hsvac")
kvac1(nv2(jspin),jspin) = lapw%k1(k,jspin)
kvac2(nv2(jspin),jspin) = lapw%k2(k,jspin)
kvac(1:2,nv2(jspin),jspin) = lapw%gvec(1:2,k,jspin)
map2(k,jspin) = nv2(jspin)
ENDDO k_loop
ENDDO
......@@ -86,20 +84,16 @@ CONTAINS
s1=MIN(SIZE(hmat,1),jspin1) !in colinear case s1=1
DO jspin2=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin
s2=MIN(SIZE(hmat,1),jspin2) !in colinear case s2=1
ipot=3
IF (jspin1==jspin2) ipot=jspin1
!---> get the wavefunctions and set up the tuuv, etc matrices
jspin=jsp
CALL vacfun(&
vacuum,DIMENSION,stars,&
jsp,input,noco,jspin1,jspin2,&
sym, cell,ivac,evac(1,1),lapw%bkpt,v%vacxy(:,:,ivac,ipot),v%vacz(:,:,:),kvac1,kvac2,nv2,&
tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
CALL vacfun(&
vacuum,stars,&
input,noco,jspin1,jspin2,&
sym, cell,ivac,evac,lapw%bkpt,v%vacxy,v%vacz,kvac,nv2,&
tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
!
!---> generate a and b coeffficients
!
IF (noco%l_noco) THEN
DO jspin = 1,input%jspins
DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
DO k = 1,lapw%nv(jspin)
gz = sign*cell%bmat(3,3)*lapw%k3(k,jspin)
i2 = map2(k,jspin)
......@@ -109,16 +103,6 @@ CONTAINS
b(k,jspin) = c_1 * CMPLX(duz(i2,jspin), gz* uz(i2,jspin) )
ENDDO
ENDDO
ELSE
DO k = 1,lapw%nv(jsp)
gz = sign*cell%bmat(3,3)*lapw%gvec(3,k,jsp)
i2 = map2(k,jsp)
th = gz*cell%z1
c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
a(k,jsp) = - c_1 * CMPLX(dudz(i2,jsp), gz*udz(i2,jsp) )
b(k,jsp) = c_1 * CMPLX(duz(i2,jsp), gz* uz(i2,jsp) )
ENDDO
ENDIF
!---> update hamiltonian and overlap matrices
IF (jspin1==jspin2) THEN
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
......@@ -153,30 +137,25 @@ CONTAINS
!Diagonal term of Overlapp matrix, Hamiltonian later
sij = CONJG(a(i,jspin2))*a(i,jspin2) + CONJG(b(i,jspin2))*b(i,jspin2)*ddnv(ik,jspin2)
IF (hmat(1,1)%l_real) THEN
smat(s1,s2)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
smat(s2,s1)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
ELSE
smat(s1,s2)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
smat(s2,s1)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
ENDIF
ENDDO
ENDIF
!---> hamiltonian update
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
DO i = mpi%n_rank+1,lapw%nv(jspin1),mpi%n_size
i0=(i-1)/mpi%n_size+1 !local column index
ik = map2(i,jspin2)
DO j = 1,MIN(i,lapw%nv(jspin1))
jk = map2(j,jspin1)
IF (jspin2>jspin1) THEN
hij = CONJG(CONJG(a(j,jspin1))* (tuuv(jk,ik)*a(i,jspin2) +tudv(jk,ik)*b(i,jspin2))&
+ CONJG(b(j,jspin1))* (tddv(jk,ik)*b(i,jspin2) +tduv(jk,ik)*a(i,jspin2)))
ELSE
hij = CONJG(a(i,jspin2))* (tuuv(ik,jk)*a(j,jspin1) +tudv(ik,jk)*b(j,jspin1))&
+ CONJG(b(i,jspin2))* (tddv(ik,jk)*b(j,jspin1) +tduv(ik,jk)*a(j,jspin1))
ENDIF
ik = map2(i,jspin1)
DO j = 1,MERGE(i,lapw%nv(jspin2),jspin1==jspin2)
jk = map2(j,jspin2)
hij = CONJG(a(i,jspin1))* (tuuv(ik,jk)*a(j,jspin2) +tudv(ik,jk)*b(j,jspin2))&
+ CONJG(b(i,jspin1))* (tddv(ik,jk)*b(j,jspin2) +tduv(ik,jk)*a(j,jspin2))
IF (hmat(1,1)%l_real) THEN
hmat(s1,s2)%data_r(j,i0) = hmat(s1,s2)%data_r(j,i0) + REAL(hij)
hmat(s2,s1)%data_r(j,i0) = hmat(s2,s1)%data_r(j,i0) + REAL(hij)
ELSE
hmat(s1,s2)%data_c(j,i0) = hmat(s1,s2)%data_c(j,i0) + hij
hmat(s2,s1)%data_c(j,i0) = hmat(s2,s1)%data_c(j,i0) + hij
ENDIF
ENDDO
ENDDO
......
This diff is collapsed.
......@@ -104,6 +104,7 @@ CONTAINS
REAL :: scale_temp, dtild_temp
REAL :: forceAllAtoms(3,atoms%nat)
CLASS(t_forcetheo),ALLOCATABLE:: forcetheo
TYPE(t_field) :: field
input=input_in
atoms_new=atoms
......@@ -207,7 +208,7 @@ CONTAINS
ALLOCATE(xmlCoreOccs(1,1,1))
CALL initWannierDefaults(wann_temp)
CALL r_inpXML(atoms_temp,obsolete_temp,vacuum_temp,input_temp,stars_temp,sliceplot_temp,&
banddos_temp,dimension_temp,forcetheo,cell_temp,sym_temp,xcpot_temp,noco_temp,&
banddos_temp,dimension_temp,forcetheo,field,cell_temp,sym_temp,xcpot_temp,noco_temp,&
oneD_temp,hybrid_temp,kpts_temp,enpara_temp,coreSpecInput_temp,wann_temp,noel_temp,&
namex_temp,relcor_temp,a1_temp,a2_temp,a3_temp,dtild_temp,xmlElectronStates,&
xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,l_kpts_temp)
......
......@@ -258,54 +258,18 @@
INTEGER, INTENT(IN) :: k1d, k2d, nvac
REAL, INTENT(IN) :: area
REAL :: zsigma, sig_b(2), tmp
INTEGER, PARAMETER :: iou = 33
INTEGER :: ios
LOGICAL :: exists, eV
CHARACTER(len=50) :: str
INQUIRE (file='apwefl',exist=exists)
IF (.NOT. exists) RETURN
OPEN (iou,file='apwefl',form='formatted',status='old')
READ (iou, '(a)') str
str = ADJUSTL (str)
! Old format: zsigma, optionally followed by
! sig_b (1) and (2)
IF (str(1:1) /= '!' .AND. str(1:1) /= '#'&
& .AND. str(1:1) /= '&' .AND. str /= '') THEN
READ (str, *) E%zsigma
READ (iou, *, iostat=ios) sig_b
IF (ios == 0) E%sig_b = sig_b
RETURN
END IF
REAL :: tmp
INTEGER :: i
! New format
ALLOCATE(E%sigEF(3*k1d, 3*k2d, nvac))
E%sigEF = 0.0
eV = .false.
REWIND(iou)
DO
READ(iou,'(a)', IOSTAT=ios) str
IF (ios /= 0) EXIT
str = ADJUSTL (str)
IF (str(1:1) == '#' .OR. str(1:1) == '!' .OR. str == '')&
& CYCLE
IF (lower_case (str(1:8)) == '&efield ') THEN
BACKSPACE (iou)
CALL read_namelist (iou, E, eV)
ELSE IF (str(1:1) == '&') THEN
CALL juDFT_error("ERROR reading 'apwefl': Unknown namelist"&
& ,calledby ="efield")
ELSE
CALL read_shape (E, str, nvac)
END IF
if (allocated(e%shapes)) then
DO i=1,SIZE(e%shapes)
CALL read_shape (E, e%shapes(i), nvac)
END DO
CLOSE (iou)
IF (eV) THEN
endif
IF (e%l_eV) THEN
E%sig_b(:) = E%sig_b/hartree_to_ev_const
E%sigEF(:,:,:) = E%sigEF/hartree_to_ev_const
END IF
......
......@@ -6,6 +6,7 @@ io/rw_symfile.f
set(fleur_F90 ${fleur_F90}
io/io_matrix.F90
io/io_hybrid.F90
io/relax_io.F90
io/eig66_da.F90
io/eig66_data.F90
io/eig66_hdf.F90
......
......@@ -14,7 +14,7 @@ MODULE m_rinpXML
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CONTAINS
SUBROUTINE r_inpXML(&
atoms,obsolete,vacuum,input,stars,sliceplot,banddos,DIMENSION,forcetheo,&
atoms,obsolete,vacuum,input,stars,sliceplot,banddos,DIMENSION,forcetheo,field,&
cell,sym,xcpot,noco,oneD,hybrid,kpts,enpara,coreSpecInput,wann,&
noel,namex,relcor,a1,a2,a3,dtild,xmlElectronStates,&
xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,&
......@@ -58,6 +58,7 @@ CONTAINS
TYPE(t_noco),INTENT(INOUT) :: noco
TYPE(t_dimension),INTENT(OUT) :: dimension
TYPE(t_enpara) ,INTENT(OUT) :: enpara
TYPE(t_field), INTENT(INOUT) :: field
CLASS(t_forcetheo),ALLOCATABLE,INTENT(OUT):: forcetheo
TYPE(t_coreSpecInput),INTENT(OUT) :: coreSpecInput
TYPE(t_wann) ,INTENT(INOUT) :: wann
......@@ -123,7 +124,7 @@ CONTAINS
REAL :: tauTemp(3,48)
REAL :: bk(3)
LOGICAL :: flipSpin, l_eV, invSym, l_qfix, relaxX, relaxY, relaxZ
LOGICAL :: l_vca, coreConfigPresent, l_enpara, l_orbcomp, tempBool, l_nocoinp
LOGICAL :: coreConfigPresent, l_enpara, l_orbcomp, tempBool, l_nocoinp
REAL :: magMom, radius, logIncrement, qsc(3), latticeScale, dr
REAL :: aTemp, zp, rmtmax, sumWeight, ldau_u(4), ldau_j(4), tempReal
REAL :: weightScale, eParamUp, eParamDown
......@@ -141,7 +142,7 @@ CONTAINS
INTEGER :: altKPointSetIndex, altKPointSetIndices(2)
LOGICAL :: ldaSpecies
REAL :: socscaleSpecies
REAL :: socscaleSpecies,b_field_mtspecies,vcaspecies
INTEGER, ALLOCATABLE :: lNumbers(:), nNumbers(:), speciesLLO(:)
INTEGER, ALLOCATABLE :: loOrderList(:)
......@@ -700,27 +701,31 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
! END IF
! Read in optional E-Field parameters
xPathA = '/fleurInput/calculationSetup/eField'
xPathA = '/fleurInput/calculationSetup/fields'
numberNodes = xmlGetNumberOfNodes(xPathA)
IF (numberNodes.EQ.1) THEN
!input%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma'))
!input%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1'))
!input%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2'))
!input%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge'))
!input%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho'))
!input%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp'))
!input%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet'))
!l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV'))
STOP 'Error: Reading input for E-Fields not yet implemented completely!'
! ALLOCATE(input%efield%sigEF(3*k1d, 3*k2d, nvac))
! input%efield%sigEF = 0.0
!IF (l_eV) THEN
! input%efield%sig_b(:) = input%efield%sig_b/hartree_to_ev_const
! input%efield%sigEF(:,:,:) = input%efield%sigEF/hartree_to_ev_const
!END IF
IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@b_field')>0) THEN
field%b_field=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'//@b_field'))
field%l_b_field=.true.
ENDIF
field%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma'))
field%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1'))
field%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2'))
field%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge'))
field%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho'))
field%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp'))
field%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet'))
field%efield%l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV'))
numberNodes=xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/shape')
ALLOCATE(field%efield%shapes(numberNodes))
DO i=1,numberNodes
WRITE(xPathB,"(a,a,i0,a)") TRIM(ADJUSTL(xpathA)),'/shape[',i,']'
field%efield%shapes(i)=TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
ENDDO
END IF
! Read in optional energy parameter limits
......@@ -1315,18 +1320,26 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
speciesNLO(iSpecies) = speciesNLO(iSpecies) + lNumCount
DEALLOCATE (lNumbers, nNumbers)
END DO
! Special switches for species
vcaspecies=0.0
WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/special'
numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
IF (numberNodes==1) THEN
vcaSpecies = evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vca_charge'))))
ENDIF
DO iType = 1, atoms%ntype
WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species'
valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
IF(TRIM(ADJUSTL(atoms%speciesName(iSpecies))).EQ.TRIM(ADJUSTL(valueString))) THEN
atoms%nz(iType) = atomicNumber
atoms%zatom(iType) = atoms%nz(iType)
IF (atoms%nz(iType).EQ.0) THEN
WRITE(*,*) 'Note: Replacing atomic number 0 by 1.0e-10 on atom type ', iType
atoms%zatom(iType) = 1.0e-10
END IF
atoms%zatom(iType)=atoms%zatom(iType)+vcaspecies
noel(iType) = namat_const(atoms%nz(iType))
atoms%zatom(iType) = atoms%nz(iType)
atoms%rmt(iType) = radius
atoms%jri(iType) = gridPoints
atoms%dx(iType) = logIncrement
......@@ -1409,6 +1422,10 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
IF (numberNodes==1) THEN
ldaSpecies = evaluateFirstBoolOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lda'))))
socscaleSpecies = evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@socscale'))))
IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@b_field_mt')>0) THEN
b_field_mtSpecies=evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@b_field_mt'))))
field%l_b_field=.TRUE.
ENDIF
ENDIF
! Explicitely provided core configurations
......@@ -1621,6 +1638,13 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
xcpot%lda_atom(iType)=ldaSpecies
END SELECT
noco%socscale(iType)=socscaleSpecies
IF (field%l_b_field) THEN
IF (.NOT.ALLOCATED(field%b_field_mt)) THEN
ALLOCATE(field%b_field_mt(atoms%ntype))
field%b_field_mt=0.0
ENDIF
field%b_field_mt(itype)=b_field_mtSpecies
ENDIF
END IF
END DO
DEALLOCATE(loOrderList)
......@@ -2119,7 +2143,7 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
!!! End of non-XML input
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CALL xmlFreeResources()
!CALL xmlFreeResources()
!WRITE(*,*) 'Reading of inp.xml file finished'
......
!--------------------------------------------------------------------------------
! Copyright (c) 2016 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_relaxio
USE m_judft
IMPLICIT NONE
PRIVATE
PUBLIC :: read_relax,write_relax
CONTAINS
SUBROUTINE write_relax(positions,forces,energies,displace)
REAL,INTENT(in):: positions(:,:,:)
REAL,INTENT(in):: forces(:,:,:)
REAL,INTENT(in):: energies(:)
REAL,INTENT(in):: displace(:,:)
INTEGER :: no_steps,n,ntype,step
No_steps=SIZE(positions,3)
ntype=SIZE(positions,2)
IF (ntype.NE.SIZE(forces,2).OR.ntype<SIZE(displace,2).OR.&
no_steps.NE.SIZE(forces,3).OR.no_steps.NE.SIZE(energies)) &
CALL judft_error("BUG in relax_io")
OPEN(765,file="relax.inp",status="replace")
WRITE(765,*) "<relaxation>"
!write current set of displacements
WRITE(765,*) " <displacements>"
DO n=1,SIZE(displace,2)
WRITE(765,"(a,i0,a,3(f15.10,1x),a)") &
' <displace na="',n,'">',displace(:,n),'</displace>'
END DO
WRITE(765,"(a)") ' </displacements>'
!Write all known old positions,forces and energies
WRITE(765,*) " <relaxation-history>"
DO step=1,no_steps
WRITE(765,"(a,f20.10,a)") ' <step energy="',energies(step),'">'
DO n=1,ntype
WRITE(765,"(a,i0,a,3(f15.10,1x),a)") &
' <pos n="',n,'">',positions(:,n,step),'</pos>'
WRITE(765,"(a,i0,a,3(f15.10,1x),a)") &
' <force n="',n,'">',forces(:,n,step),'</force>'
END DO
WRITE(765,"(a)") ' </step>'
ENDDO
WRITE(765,*) " </relaxation-history>"
WRITE(765,*) "</relaxation>"
CLOSE(765)
END SUBROUTINE write_relax
SUBROUTINE read_relax(positions,forces,energies)
USE m_xmlIntWrapFort
USE m_calculator
REAL,INTENT(OUT),ALLOCATABLE:: positions(:,:,:)
REAL,INTENT(OUT),ALLOCATABLE:: forces(:,:,:)
REAL,INTENT(OUT),ALLOCATABLE:: energies(:)
INTEGER:: no_steps
INTEGER:: ntype,step,n
CHARACTER(len=50):: path,p,str
no_steps=xmlGetNumberOfNodes('/fleurInput/relaxation/relaxation-history/step')
IF (no_steps==0) THEN
ALLOCATE(positions(0,0,0),forces(0,0,0),energies(0))
RETURN
ELSE
ntype=xmlGetNumberOfNodes('/fleurInput/relaxation/relaxation-history/step[1]/pos')
END IF
ALLOCATE(positions(3,ntype,no_steps))
ALLOCATE(forces(3,ntype,no_steps))
ALLOCATE(energies(no_steps))
DO step=1,no_steps
WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/relaxation-history/step[',step,']'
energies(step)=evaluateFirstOnly(xmlGetAttributeValue(TRIM(path)//"/energy"))
DO n=1,ntype
WRITE(p,"(a,a,i0,a)") TRIM(path),"/pos[",n,"]"
str=xmlGetAttributeValue(p)
positions(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
WRITE(p,"(a,a,i0,a)") TRIM(path),"/force[",n,"]"
str=xmlGetAttributeValue(p)
forces(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
ENDDO
END DO
END SUBROUTINE read_relax
SUBROUTINE read_displacements(cell,atoms)
USE m_xmlIntWrapFort
USE m_calculator
USE m_types
TYPE(t_cell),INTENT(INOUT) :: cell
TYPE(t_atoms),INTENT(INOUT):: atoms
INTEGER:: n,na
REAL :: d(3)
CHARACTER(len=50):: path,str
IF (xmlGetNumberOfNodes('/fleurInput/relaxation/displacements').NE.0) THEN
!read displacements and apply to positions
DO n=1,xmlGetNumberOfNodes('/fleurInput/relaxation/displacements/displace')
WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/displacements/displace[',n,']'
str=xmlGetAttributeValue(path)
d=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
na=INT(evaluateFirstOnly(xmlGetAttributeValue(TRIM(path)//"/@na")))
IF (na>SIZE(atoms%taual,2)) CALL judft_error&
("Wrong number of displacements. na can not be larger than number of atoms in setup")
atoms%taual(:,na)=atoms%taual(:,na)+d
atoms%pos(:,na) = matmul(cell%amat,atoms%taual(:,na))
END DO
END IF
END SUBROUTINE read_displacements
END MODULE m_relaxio
......@@ -52,7 +52,7 @@
<xsd:element maxOccurs="1" minOccurs="0" name="geometryOptimization" type="GeometryOptimizerType"/>
<xsd:element maxOccurs="1" minOccurs="0" name="ldaU" type="LdaUGeneralType"/>
<xsd:element maxOccurs="1" minOccurs="0" name="spinSpiralQPointMesh" type="SpinSpiralQPointMeshType"/>
<xsd:element maxOccurs="1" minOccurs="0" name="eField" type="EFieldType"/>
<xsd:element maxOccurs="1" minOccurs="0" name="fields" type="FieldType"/>
<xsd:element maxOccurs="1" minOccurs="0" name="energyParameterLimits" type="EnergyParameterLimitsType"/>
</xsd:all>
</xsd:complexType>
......@@ -62,7 +62,11 @@
<xsd:attribute name="elup" type="xsd:string"/>
</xsd:complexType>
<xsd:complexType name="EFieldType">
<xsd:complexType name="FieldType">
<xsd:sequence>
<xsd:element maxOccurs="unbounded" minOccurs="0" name="shape" type="xsd:string" />
</xsd:sequence>
<xsd:attribute default="0.0" name="b_field" type="xsd:string" use="optional"/>
<xsd:attribute default="10.0" name="zsigma" type="xsd:string" use="optional"/>
<xsd:attribute default="0.0" name="sig_b_1" type="xsd:string" use="optional"/>
<xsd:attribute default="0.0" name="sig_b_2" type="xsd:string" use="optional"/>
......@@ -356,6 +360,8 @@
<xsd:complexType name="SpecialType">
<xsd:attribute name="lda" type="FleurBool" use="optional" default="F"/>
<xsd:attribute name="socscale" type="ZeroToOneNumberType" use="optional" default="1.0"/>
<xsd:attribute name="b_field_mt" type="xsd:string" use="optional" default="0.0"/>
<xsd:attribute name="vca_charge" type="xsd:string" use="optional" default="0.0"/>
</xsd:complexType>
......
......@@ -172,7 +172,7 @@
a2 = 0.0
a3 = 0.0
CALL r_inpXML(&
atoms,obsolete,vacuum,input,stars,sliceplot,banddos,DIMENSION,forcetheo,&
atoms,obsolete,vacuum,input,stars,sliceplot,banddos,DIMENSION,forcetheo,field,&
cell,sym,xcpot,noco,oneD,hybrid,kpts,enpara,coreSpecInput,wann,&
noel,namex,relcor,a1,a2,a3,dtild,xmlElectronStates,&
xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,&
......
......@@ -24,6 +24,8 @@ MODULE m_types_field
LOGICAL :: autocomp = .TRUE. ! Auto-compensate film charge
LOGICAL :: dirichlet = .FALSE. ! Dirichlet vs. Neumann boundary cond.
LOGICAL :: l_dirichlet_coeff = .FALSE. ! For MPI, true if C1/C2 set
LOGICAL :: l_eV =.FALSE. !Input in eV
CHARACTER(len=50),ALLOCATABLE :: shapes(:)
END TYPE t_efield
TYPE t_field
......
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