Commit 7588283e authored by Daniel Wortmann's avatar Daniel Wortmann

Cleanup of forces part, new input, etc...

parent 3a758e0e
set(fleur_F77 ${fleur_F77}
force/relax.F
#force/relax.F
)
set(fleur_F90 ${fleur_F90}
force/bfgs0.f90
force/bfgs.f90
#force/bfgs0.f90
#force/bfgs.f90
force/force_b8.f90
force/force_0.f90
#force/force_0.f90
force/force_a12.f90
force/force_a21.F90
force/force_a21_U.f90
......@@ -17,7 +17,7 @@ force/force_a8.F90
force/force_b8.f90
force/force_sf.F90
force/force_w.F90
force/geo.f90
#force/geo.f90
force/stern.f90
force/relaxation.F90
)
MODULE m_forcew
! ************************************************************
! Printing force components
! ************************************************************
CONTAINS
!--------------------------------------------------------------------------------
! 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_forcew
! ************************************************************
! Printing force components
! ************************************************************
CONTAINS
SUBROUTINE force_w(mpi,input,atoms,sym,results,cell,oneD,vacuum)
USE m_geo
USE m_relax
USE m_types
USE m_xmlOutput
use m_relaxation
USE m_relaxation
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_results),INTENT(IN) :: results
TYPE(t_results),INTENT(INOUT):: results
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_input),INTENT(IN) :: input
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_vacuum),INTENT(IN) :: vacuum
! ..
! .. Local Scalars ..
! ..
! .. Local Scalars ..
REAL,PARAMETER:: zero=0.0
REAL sum
INTEGER i,jsp,n,nat1,ierr
REAL eps_force
LOGICAL :: l_new,l_relax
! ..
! .. Local Arrays ..
! ..
! .. Local Arrays ..
REAL forcetot(3,atoms%ntype)
CHARACTER(LEN=20) :: attributes(7)
!
! write spin-dependent forces
!
!
! write spin-dependent forces
!
IF (mpi%irank==0) THEN
nat1 = 1
DO n = 1,atoms%ntype
......@@ -42,17 +45,17 @@
& (results%force(i,n,jsp),i=1,3)
END DO
END IF
8000 FORMAT ('SPIN-',i1,1x,'FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,&
8000 FORMAT ('SPIN-',i1,1x,'FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,&
& 3x,'Y=',f7.3,3x,'Z=',f7.3,5x,' FX_SP =',f9.6,' FY_SP =',&
& f9.6,' FZ_SP =',f9.6)
ENDIF
nat1 = nat1 + atoms%neq(n)
END DO
!
! write total forces
!
!
! write total forces
!
WRITE (6,8005)
8005 FORMAT (/,' ***** TOTAL FORCES ON ATOMS ***** ',/)
8005 FORMAT (/,' ***** TOTAL FORCES ON ATOMS ***** ',/)
IF (input%l_f) CALL openXMLElement('totalForcesOnRepresentativeAtoms',(/'units'/),(/'Htr/bohr'/))
nat1 = 1
DO n = 1,atoms%ntype
......@@ -68,7 +71,7 @@
WRITE (6,FMT=8010) n, (atoms%pos(i,nat1),i=1,3),&
& (forcetot(i,n),i=1,3)
8010 FORMAT (' TOTAL FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,3x,'Y=',&
8010 FORMAT (' TOTAL FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,3x,'Y=',&
& f7.3,3x,'Z=',f7.3,/,22x,' FX_TOT=',f9.6,&
& ' FY_TOT=',f9.6,' FZ_TOT=',f9.6)
......@@ -82,38 +85,28 @@
IF (input%l_f) THEN
CALL writeXMLElementFormPoly('forceTotal',(/'atomType','x ','y ','z ',&
'F_x ','F_y ','F_z '/),attributes,&
reshape((/8,1,1,1,3,3,3,6,12,12,12,12,12,12/),(/7,2/)))
RESHAPE((/8,1,1,1,3,3,3,6,12,12,12,12,12,12/),(/7,2/)))
END IF
END IF
nat1 = nat1 + atoms%neq(n)
END DO
IF (input%l_f) CALL closeXMLElement('totalForcesOnRepresentativeAtoms')
sum=0.0
DO n = 1,atoms%ntype
IF (atoms%l_geo(n)) THEN
DO i = 1,3
sum = max(sum,(forcetot(i,n) - results%force_old(i,n))**2)
ENDDO
ENDIF
ENDDO
sum=sqrt(sum)
!-roa
eps_force=0.00001
open(88,file='eps_force',form='formatted',status='old',err=188)
read(88,'(f20.8)') eps_force
close (88)
188 continue
!-roa
WRITE (6,8020) eps_force,sum
8020 FORMAT ('eps_force=',f8.5,'max=',f8.5)
!Check convergence of force by comparing force with old_force
sum=MAXVAL(ABS(forcetot - results%force_old))
results%force_old(:,:)=forcetot !Store for next iteration
results%force=0.0
l_relax=sum<input%force_converged
IF (.NOT.l_relax) THEN
WRITE (6,8020) input%force_converged,sum
8020 FORMAT ('No new postions, force convergence required=',f8.5,'max force distance=',f8.5)
END IF
ENDIF
l_relax=sum<eps_force
#ifdef CPP_MPI
CALL MPI_BCAST(l_relax,1,MPI_LOGICAL,0,ierr)
#endif
IF (l_relax.and.input%l_f) CALL relaxation(mpi,input,atoms,cell,sym,forcetot,results%tote)
IF (l_relax.AND.input%l_f) CALL relaxation(mpi,input,atoms,cell,sym,forcetot,results%tote)
END SUBROUTINE force_w
END MODULE m_forcew
!--------------------------------------------------------------------------------
! 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_relaxation
USE m_judft
IMPLICIT NONE
PRIVATE
integer:: input_force_relax=3
public relaxation
PUBLIC relaxation !This is the interface. Below there are internal subroutines for bfgs, simple mixing, CG ...
CONTAINS
SUBROUTINE relaxation(mpi,input,atoms,cell,sym,force_new,energies_new)
!This routine uses the current force,energies and atomic positions to
!generate a displacement in a relaxation step.
!The history is taken into account by read_relax from m_relaxio
!After generating new positions the code stops
USE m_types
use m_relaxio
use m_broyd_io
USE m_relaxio
USE m_broyd_io
#ifdef CPP_MPI
INCLUDE 'mpif.h'
#endif
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_input),INTENT(IN) :: input
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
REAL,INTENT(in) :: force_new(:,:),energies_new
REAL,INTENT(in) :: force_new(:,:),energies_new !data for this iteration
REAL,ALLOCATABLE :: pos(:,:,:),force(:,:,:),energies(:)
REAL,ALLOCATABLE :: displace(:,:),old_displace(:,:)
REAL :: alpha
INTEGER :: n
INTEGER :: n,ierr
LOGICAL :: l_conv
IF (mpi%irank==0) THEN
ALLOCATE(pos(3,atoms%ntype,1));
......@@ -36,32 +47,44 @@ CONTAINS
CALL read_relax(pos,force,energies)
!determine new positions
IF (SIZE(energies)==1.OR.input_force_relax==0) THEN
IF (SIZE(energies)==1.OR.input%forcemix==0) THEN
!no history present simple step
! choose a reasonable first guess for scaling
! this choice is based on a Debye temperature of 330K;
! modify as needed
alpha = (250.0/(MAXVAL(atoms%zatom)*input%xa))*((330./input%thetad)**2)
CALL simple_step(alpha,force,displace)
ELSEIF (input_force_relax==1) THEN
!alpha = (250.0/(MAXVAL(atoms%zatom)*input%xa))*((330./input%thetad)**2)
CALL simple_step(input%forcealpha,force,displace)
ELSEIF (input%forcemix==1) THEN
CALL simple_cg(pos,force,displace)
ELSE
CALL simple_bfgs(pos,force,displace)
ENDIF
CALL read_displacements(atoms,old_displace)
!Check for convergence of forces/displacements
l_conv=.TRUE.
DO n=1,atoms%ntype
PRINT *,"OD:",old_displace(:,n)
PRINT *,"ND:",displace(:,n)
END DO
IF (DOT_PRODUCT(force(:,n,SIZE(force,3)),force(:,n,SIZE(force,3)))>input%epsforce**2) l_conv=.FALSE.
IF (DOT_PRODUCT(displace(:,n),displace(:,n))>input%epsforce**2) l_conv=.FALSE.
ENDDO
!New displacements relative to positions in inp.xml
CALL read_displacements(atoms,old_displace)
displace=displace+old_displace
!Write file
CALL write_relax(pos,force,energies,displace)
ENDIF
#ifdef CPP_MPI
CALL MPI_BCAST(l_conv,1,MPI_LOGICAL,0,ierr)
#endif
IF (l_conv) THEN
CALL judft_end("Structual relaxation: Done",0)
ELSE
CALL resetBroydenHistory()
CALL judft_end("Structual relaxation done",0)
CALL judft_end("Structual relaxation: new displacements generated",0)
END IF
END SUBROUTINE relaxation
......@@ -83,16 +106,16 @@ CONTAINS
!-----------------------------------------------
IMPLICIT NONE
REAL,INTENT(in) :: pos(:,:,:),force(:,:,:)
real,INTENT(OUT) :: shift(:,:)
REAL,INTENT(OUT) :: shift(:,:)
INTEGER :: n,i,j,hist_length,n_force
REAL,ALLOCATABLE:: h(:,:)
REAL,ALLOCATABLE:: p(:),y(:),v(:)
REAL :: py,yy,gamma
n_force=3*size(pos,2)
allocate(h(n_force,n_force))
allocate(p(n_force),y(n_force),v(n_force))
n_force=3*SIZE(pos,2)
ALLOCATE(h(n_force,n_force))
ALLOCATE(p(n_force),y(n_force),v(n_force))
!calculate approx. Hessian
!initialize H
......@@ -101,15 +124,15 @@ CONTAINS
h(n,n) = 1.0
ENDDO
!loop over all iterations (including current)
hist_length=size(pos,3)
hist_length=SIZE(pos,3)
DO n = 2,hist_length
! differences
p(:) = RESHAPE(pos(:,:,n)-pos(:,:,n-1),(/SIZE(p)/))
y(:) = RESHAPE(force(:,:,n)-force(:,:,n-1),(/SIZE(p)/))
! get necessary inner products and H|y>
py = dot_PRODUCT(p,y)
py = DOT_PRODUCT(p,y)
v = MATMUL(y,h)
yy = dot_PRODUCT(y,v)
yy = DOT_PRODUCT(y,v)
!check that update will leave h positive definite;
IF (py <= 0.0) THEN
WRITE (6,*) ' bfgs: <p|y> < 0'
......@@ -135,14 +158,14 @@ CONTAINS
ENDIF
ENDDO
y(:) = RESHAPE(force(:,:,hist_length),(/SIZE(p)/))
shift = reshape(MATMUL(y,h),shape(shift))
shift = RESHAPE(MATMUL(y,h),SHAPE(shift))
END SUBROUTINE simple_bfgs
SUBROUTINE simple_cg(pos,force,shift)
!-----------------------------------------------
IMPLICIT NONE
REAL,intent(in) :: pos(:,:,:),force(:,:,:)
real,INTENT(OUT) :: shift(:,:)
REAL,INTENT(in) :: pos(:,:,:),force(:,:,:)
REAL,INTENT(OUT) :: shift(:,:)
REAL :: f1(3,SIZE(pos,2)),f2(3,SIZE(pos,2))
INTEGER :: n_old
......
......@@ -155,7 +155,7 @@
banddos%ndir = 0 ; vacuum%layers = 0 ; atoms%nflip(:) = 1 ; vacuum%izlay(:,:) = 0
banddos%e_mcd_lo = -10.0 ; banddos%e_mcd_up = 0.0
atoms%lda_u%l = -1 ; atoms%relax(1:2,:) = 1 ; atoms%relax(:,:) = 1
input%epsdisp = 0.00001 ; input%epsforce = 0.00001 ; input%xa = 2.0 ; input%thetad = 330.0
input%epsdisp = 0.00001 ; input%epsforce = 0.00001 ; input%forcealpha = 1.0
sliceplot%e1s = 0.0 ; sliceplot%e2s = 0.0 ; banddos%e1_dos = 0.5 ; banddos%e2_dos = -0.5 ; input%tkb = 0.001
banddos%sig_dos = 0.015 ; vacuum%tworkf = 0.0 ; input%scaleCell = 1.0 ; scpos = 1.0
input%scaleA1 = 1.0 ; input%scaleA2 = 1.0 ; input%scaleC = 1.0
......
......@@ -666,11 +666,11 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
IF (numberNodes.EQ.1) THEN
input%l_f = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_f'))
input%xa = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@xa'))
input%thetad = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@thetad'))
input%forcealpha = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@forcealpha'))
input%epsdisp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsdisp'))
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
......
......@@ -5,6 +5,9 @@
!--------------------------------------------------------------------------------
MODULE m_relaxio
!This module handles IO to the relax.xml file
!The writing is done directly to relax.xml
!The reading uses the libxml interface to inp.xml. Hence the relax.xml has to be included here.
USE m_judft
IMPLICIT NONE
PRIVATE
......@@ -19,7 +22,7 @@ CONTAINS
INTEGER :: no_steps,n,ntype,step
No_steps=SIZE(positions,3)
ntype=SIZE(positions,2)
IF (ntype.NE.SIZE(forces,2).OR.ntype.ne.SIZE(displace,2).OR.&
IF (ntype.NE.SIZE(forces,2).OR.ntype.NE.SIZE(displace,2).OR.&
no_steps.NE.SIZE(forces,3).OR.no_steps.NE.SIZE(energies))THEN
CALL judft_error("BUG in relax_io")
ENDIF
......@@ -68,15 +71,15 @@ CONTAINS
IF (ALLOCATED(positions)) THEN
!Assume that we got already a set of positions, forces, energy and extend that list
rtmp=positions
deallocate(positions)
DEALLOCATE(positions)
ALLOCATE(positions(3,ntype,no_steps+SIZE(rtmp,3)))
positions(:,:,no_steps+1:)=rtmp
rtmp=forces
deallocate(forces)
DEALLOCATE(forces)
ALLOCATE(forces(3,ntype,no_steps+SIZE(rtmp,3)))
forces(:,:,no_steps+1:)=rtmp
rtmp(1,1,:)=energies
deallocate(energies)
DEALLOCATE(energies)
ALLOCATE(energies(no_steps+SIZE(rtmp,3)))
energies(no_steps+1:)=rtmp(1,1,:)
ELSE
......@@ -104,7 +107,7 @@ CONTAINS
TYPE(t_atoms),INTENT(in)::atoms
REAL,INTENT(out)::disp(:,:)
CHARACTER(len=50):: path,str
integer :: n
INTEGER :: n
disp=0.0
IF (xmlGetNumberOfNodes('/fleurInput/relaxation/displacements')==0) RETURN
!read displacements and apply to positions
......@@ -119,7 +122,7 @@ CONTAINS
SUBROUTINE apply_displacements(cell,input,vacuum,oneD,sym,noco,atoms)
USE m_types
USE m_chkmt
use m_mapatom
USE m_mapatom
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_cell),INTENT(IN) :: cell
......
......@@ -220,8 +220,8 @@ SUBROUTINE w_inpXML(&
WRITE (fileNum,180) input%gw,input%secvar
! <geometryOptimization l_f="F" xa="2.00000" thetad="330.00000" epsdisp="0.00001" epsforce="0.00001"/>
190 FORMAT(' <geometryOptimization l_f="',l1,'" xa="',f0.8,'" thetad="',f0.8,'" epsdisp="',f0.8,'" epsforce="',f0.8,'"/>')
WRITE (fileNum,190) input%l_f,input%xa,input%thetad,input%epsdisp,input%epsforce
190 FORMAT(' <geometryOptimization l_f="',l1,'" forcealpha="',f0.8,'" forcemix="',i0,'" epsdisp="',f0.8,'" epsforce="',f0.8,'"/>')
WRITE (fileNum,190) input%l_f,input%forcealpha,input%forcemix,input%epsdisp,input%epsforce
IF(input%gauss.AND.input%tria) THEN
STOP 'Error: bz integration modes gauss AND tria selected!'
......
......@@ -587,10 +587,10 @@
<xsd:complexType name="GeometryOptimizerType">
<xsd:attribute name="l_f" type="FleurBool" use="required"/>
<xsd:attribute name="xa" type="xsd:string" use="required"/>
<xsd:attribute name="thetad" type="xsd:string" use="required"/>
<xsd:attribute name="epsdisp" type="xsd:string" use="required"/>
<xsd:attribute name="epsforce" type="xsd:string" use="required"/>
<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:complexType>
......
......@@ -49,7 +49,6 @@ CONTAINS
USE m_eigen
USE m_eigenso
USE m_fermie
USE m_force0
USE m_cdngen
USE m_totale
USE m_potdis
......@@ -330,7 +329,6 @@ CONTAINS
CYCLE forcetheoloop
ENDIF
CALL force_0(results)! ----> initialise force_old
!!$ !+Wannier functions
!!$ IF ((input%l_wann).AND.(.NOT.wann%l_bs_comf)) THEN
......
......@@ -152,8 +152,8 @@ CONTAINS
! neigd2 = dimension%neigd
IF (noco%l_soc.AND.(.NOT.noco%l_noco)) neigd2 = 2*neigd2
ALLOCATE (thisResults%force(3,atoms%ntype,input%jspins))
ALLOCATE (thisResults%force_old(3,atoms%ntype))
ALLOCATE (thisResults%force(3,atoms%ntype,input%jspins));thisResults%force=0.0
ALLOCATE (thisResults%force_old(3,atoms%ntype));thisResults%force_old=0.0
ALLOCATE (thisResults%w_iks(neigd2,kpts%nkpt,input%jspins))
ALLOCATE (thisResults%neig(kpts%nkpt,input%jspins))
ALLOCATE (thisResults%eig(neigd2,kpts%nkpt,input%jspins))
......
......@@ -374,10 +374,11 @@ MODULE m_types_setup
INTEGER :: gw
INTEGER :: gw_neigd
INTEGER :: qfix
REAL :: xa !< mixing parameter for geometry optimzer
REAL :: thetad !< Debey temperature for first step of geometry optimzer
REAL :: forcealpha !< mixing parameter for geometry optimzer
REAL :: epsdisp !< minimal displacement. If all displacements are < epsdisp stop
REAL :: epsforce !< minimal force. If all forces <epsforce stop
REAL :: force_converged=0.00001
REAL :: forcemix=3
REAL :: delgau
REAL :: alpha
REAL :: preconditioning_param
......
......@@ -45,7 +45,7 @@ module m_VYukawaFilm
! PSEUDO-CHARGE DENSITY
call psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, psq )
call psqpw( mpi, atoms, sphhar, stars, vacuum, cell, input, sym, oneD, den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, psq )
! VACUUM POTENTIAL
......
......@@ -70,7 +70,7 @@ CONTAINS
atoms_tmp%zatom=0.0
CALL psqpw(mpi,&
atoms_tmp,sphhar,stars,vacuum,&
DIMENSION,cell,input,sym,oneD,&
cell,input,sym,oneD,&
qpw,rho,(/0.,0./),.TRUE.,2,psq)
!put pseudo charge on real-space grid
......
......@@ -16,7 +16,7 @@ module m_psqpw
contains
subroutine psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, &
subroutine psqpw( mpi, atoms, sphhar, stars, vacuum, cell, input, sym, oneD, &
& qpw, rho, rht, l_xyav, potdenType, psq )
#include"cpp_double.h"
......@@ -37,7 +37,6 @@ contains
type(t_sphhar), intent(in) :: sphhar
type(t_stars), intent(in) :: stars
type(t_vacuum), intent(in) :: vacuum
type(t_dimension), intent(in) :: dimension
type(t_cell), intent(in) :: cell
type(t_input), intent(in) :: input
type(t_sym), intent(in) :: sym
......@@ -56,7 +55,7 @@ contains
complex :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
real :: q2(vacuum%nmzd)
real :: pn(0:atoms%lmaxd,atoms%ntype)
real :: aj(0:atoms%lmaxd+DIMENSION%ncvd+1)
real :: aj(0:atoms%lmaxd+maxval(atoms%ncv)+1)
real :: rht1(vacuum%nmz)
real, allocatable, dimension(:) :: il, kl
real :: g0(atoms%ntype)
......
......@@ -75,7 +75,7 @@ contains
! PSEUDO-CHARGE DENSITY COEFFICIENTS
call timestart( "psqpw" )
call psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, &
call psqpw( mpi, atoms, sphhar, stars, vacuum, cell, input, sym, oneD, &
den%pw(:,ispin), den%mt(:,:,:,ispin), den%vacz(:,:,ispin), .false., vCoul%potdenType, psq )
call timestop( "psqpw" )
......
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