Redesign of relaxation part

force/CMakeLists.txt

@@ -16,7 +16,8 @@ force/force_a4_add.f90
 force/force_a8.F90
 force/force_b8.f90
 force/force_sf.F90
-force/force_w.f90
+force/force_w.F90
 force/geo.f90
 force/stern.f90
+force/relaxation.F90
 )

force/force_w.f90 → force/force_w.F90

@@ -3,14 +3,14 @@
 ! Printing force components
 ! ************************************************************
       CONTAINS
-      SUBROUTINE force_w(&
-     &                   input,atoms,sym,results,cell,oneD,vacuum)
+        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
       IMPLICIT NONE
-
+      TYPE(t_mpi),INTENT(IN)       :: mpi
       TYPE(t_results),INTENT(IN)   :: results
       TYPE(t_oneD),INTENT(IN)      :: oneD
       TYPE(t_input),INTENT(IN)     :: input
@@ -24,7 +24,7 @@
       REAL sum
       INTEGER i,jsp,n,nat1,ierr
       REAL eps_force
-      LOGICAL :: l_new
+      LOGICAL :: l_new,l_relax
 !     ..
 !     .. Local Arrays ..
       REAL forcetot(3,atoms%ntype)
@@ -32,6 +32,7 @@
 !
 !     write spin-dependent forces
 !
+      IF (mpi%irank==0) THEN
       nat1 = 1
       DO n = 1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
@@ -107,19 +108,12 @@
  
       WRITE (6,8020) eps_force,sum
  8020 FORMAT ('eps_force=',f8.5,'max=',f8.5)
-
-      INQUIRE(file ="relax_inp",exist= l_new)
-      IF (l_new) THEN
-        CALL relax(input%film,atoms%pos,atoms%neq,sym%mrot,sym%tau,cell%amat,cell%bmat,atoms%ngopr,sym%invtab&
-     &        ,forcetot)
-      ELSE
-
-         IF ((sum<eps_force).AND.input%l_f) THEN
-            CALL geo(&
-     &           atoms,sym,cell,&
-     &           oneD,vacuum,input,results%tote,forcetot)
-         END IF
-      ENDIF
-
-      END SUBROUTINE force_w
-      END MODULE m_forcew
+   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)
+      
+ END SUBROUTINE force_w
+END MODULE m_forcew

force/relaxation.F90

+MODULE m_relaxation
+  USE m_judft
+  IMPLICIT NONE
+  PRIVATE
+  integer:: input_force_relax=3
+  public relaxation
+  
+CONTAINS
+  
+  SUBROUTINE relaxation(mpi,input,atoms,cell,sym,force_new,energies_new)
+    USE m_types
+    use m_relaxio
+    use m_broyd_io
+    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,ALLOCATABLE :: pos(:,:,:),force(:,:,:),energies(:)
+    REAL,ALLOCATABLE :: displace(:,:),old_displace(:,:)
+    REAL             :: alpha
+    INTEGER          :: n
+  
+    IF (mpi%irank==0) THEN
+       ALLOCATE(pos(3,atoms%ntype,1)); 
+       DO n=1,atoms%ntype
+          pos(:,n,1)=atoms%taual(:,SUM(atoms%neq(:n-1))+1)
+       END DO
+       ALLOCATE(force(3,atoms%ntype,1)); force(:,:,1)=force_new
+       ALLOCATE(energies(1));energies(1)=energies_new
+       ALLOCATE(displace(3,atoms%ntype),old_displace(3,atoms%ntype))
+    
+    ! add history 
+    CALL read_relax(pos,force,energies)
+    
+    !determine new positions
+    IF (SIZE(energies)==1.OR.input_force_relax==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
+       CALL simple_cg(pos,force,displace)
+    ELSE
+       CALL simple_bfgs(pos,force,displace)
+    ENDIF
+    
+    CALL read_displacements(atoms,old_displace)
+    DO n=1,atoms%ntype
+       PRINT *,"OD:",old_displace(:,n)
+       PRINT *,"ND:",displace(:,n)
+    END DO
+
+    displace=displace+old_displace
+    
+    !Write file
+    CALL write_relax(pos,force,energies,displace)
+ ENDIF
+ CALL resetBroydenHistory()
+ CALL judft_end("Structual relaxation done",0)
+  END SUBROUTINE relaxation
+  
+
+
+  SUBROUTINE simple_step(alpha,force,displace)
+    !-----------------------------------------------
+    IMPLICIT NONE
+    REAL,INTENT(in)  :: alpha
+    REAL,INTENT(in)  :: force(:,:,:)
+    REAL,INTENT(OUT) :: displace(:,:)
+    
+    
+    displace = alpha*force(:,:,SIZE(force,3))
+  END SUBROUTINE simple_step
+  
+  SUBROUTINE simple_bfgs(pos,force,shift)
+    !-----------------------------------------------
+    !  Simple BFGS method to calculate shift out of old positions and forces
+    !-----------------------------------------------
+    IMPLICIT NONE
+    REAL,INTENT(in)  :: pos(:,:,:),force(:,:,:)
+    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))
+
+    !calculate approx. Hessian
+    !initialize H
+    h = 0.0
+    DO n = 1,n_force
+       h(n,n) = 1.0
+    ENDDO
+    !loop over all iterations (including current)
+    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)
+       v = MATMUL(y,h)
+       yy = dot_PRODUCT(y,v)
+       !check that update will leave h positive definite;
+       IF (py <= 0.0) THEN
+          WRITE (6,*) '  bfgs: <p|y> < 0'
+          WRITE (6,*) '  check convergence of forces'
+          !Starting over with initial hessian
+          h = 0.0
+          DO j = 1,n_force
+             h(j,j) = 1.0
+          ENDDO
+          CYCLE 
+       ELSE
+          !update h
+          IF (n == 2) THEN
+             gamma = py/yy
+          ELSE
+             gamma = 1.0
+          ENDIF
+          DO j = 1,n_force
+             DO i = 1,n_force
+                h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma + (1.+gamma*yy/py)*p(i)*p(j)/py
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDDO
+    y(:) = RESHAPE(force(:,:,hist_length),(/SIZE(p)/))
+    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                :: f1(3,SIZE(pos,2)),f2(3,SIZE(pos,2))
+    INTEGER             :: n_old
+    
+    n_old = SIZE(pos,3)-1
+    
+    f1 = (force(:,:,n_old+1)-force(:,:,n_old))/(pos(:,:,n_old+1)-pos(:,:,n_old))
+    f2 = force(:,:,n_old+1)-f1*pos(:,:,n_old+1)
+    shift = -1.*f2/f1-force(:,:,n_old+1)
+  END SUBROUTINE simple_cg
+END MODULE m_relaxation

global/chkmt.f90

@@ -15,10 +15,9 @@
 !                         GM'16
 !---------------------------------------------------------------------
       CONTAINS
-      SUBROUTINE chkmt(&
-     &                 atoms,input,vacuum,cell,oneD,&
-     &                 l_gga,noel,l_test,&
-     &                 kmax,dtild,dvac1,lmax1,jri1,rmt1,dx1)
+        SUBROUTINE chkmt(atoms,input,vacuum,cell,oneD,l_test,&
+             l_gga,noel, kmax,dtild,dvac1,lmax1,jri1,rmt1,dx1,&!optional, if l_gga and ... are present suggestions are calculated
+             overlap)!this is optional, if present and l_test the routine returns the overlaps and does not stop
 
       USE m_types
       USE m_sort
@@ -32,13 +31,15 @@
       TYPE(t_vacuum),INTENT(IN):: vacuum
       TYPE(t_cell),INTENT(IN)  :: cell
       TYPE(t_oneD),INTENT(IN)  :: oneD
-      CHARACTER*3, INTENT (IN) :: noel(atoms%ntype)
-      LOGICAL, INTENT (IN)     :: l_gga,l_test
-      REAL,    INTENT (OUT)    :: kmax,dtild,dvac1
+      CHARACTER*3, INTENT (IN),OPTIONAL :: noel(atoms%ntype)
+      LOGICAL, INTENT (IN),OPTIONAL     :: l_gga
+      LOGICAL, INTENT (IN)     ::l_test
+      REAL,    INTENT (OUT),OPTIONAL    :: kmax,dtild,dvac1
 !     ..
 !     .. Array Arguments ..
-      INTEGER, INTENT (OUT)    :: lmax1(atoms%ntype),jri1(atoms%ntype)
-      REAL,    INTENT (OUT)    :: rmt1(atoms%ntype),dx1(atoms%ntype)
+      INTEGER, INTENT (OUT),OPTIONAL    :: lmax1(atoms%ntype),jri1(atoms%ntype)
+      REAL,    INTENT (OUT),OPTIONAL    :: rmt1(atoms%ntype),dx1(atoms%ntype)
+      REAL,OPTIONAL,INTENT(OUT):: overlap(0:atoms%ntype,atoms%ntype)
 !     ..
 !     .. Local Scalars ..
       INTEGER na,n
@@ -92,6 +93,7 @@
 
 
 !     0. Do initializations and set some constants
+      if (present(overlap)) overlap=0.0
 
       rmtMaxDefault = 2.8
       rmtMax = rmtMaxDefault
@@ -291,6 +293,7 @@
          END IF
       END DO
 
+      IF (PRESENT(l_gga)) THEN
 !        Sort distances and set MT radii for the atoms
 
       CALL sort(sortedDistList,nearestAtomDists)
@@ -405,6 +408,7 @@
       WRITE (6,'("k_max =",f8.5)') rkm
       WRITE (6,'("G_max =",f8.5)') 3*rkm
       kmax = rkm
+      ENDIF
 
 !     7. Test old MT radii
 
@@ -415,6 +419,7 @@
                k = nearestNeighbors(j,i)
                IF (atoms%rmt(i)+atoms%rmt(k).GE.nearestNeighborDists(j,i)) THEN
                   error = .TRUE.
+                  IF (PRESENT(overlap)) overlap(i,k)=atoms%rmt(i)+atoms%rmt(k)-nearestNeighborDists(j,i)
                   WRITE(6,240) i,k,nearestNeighborDists(j,i),atoms%rmt(i),atoms%rmt(k)
                END IF
             END DO
@@ -434,13 +439,14 @@
                          ((atoms%pos(3,iAtom)-atoms%rmt(i)).LT.-vacuum%dvac/2.)) THEN
                         error=.TRUE.
                         WRITE(6,241) i ,na
+                        IF (PRESENT(overlap)) overlap(0,i)=MAX(atoms%pos(3,iAtom)+atoms%rmt(i)-vacuum%dvac/2.,atoms%pos(3,iAtom)-atoms%rmt(i)+vacuum%dvac/2.)
                         WRITE(6,*) atoms%pos(3,iAtom),atoms%rmt(i),vacuum%dvac/2.
                      ENDIF
                   ENDIF
                END DO
             END IF
          END DO
-         IF (error) CALL juDFT_error("Error checking M.T. radii",calledby ="chkmt")
+         IF (error.AND..NOT.PRESENT(overlap)) CALL juDFT_error("Error checking M.T. radii",calledby ="chkmt")
       END IF
 
       DEALLOCATE(nearestNeighbors,numNearestNeighbors,nearestNeighborDists)

init/old_inp/inped.F90

@@ -349,7 +349,7 @@
           !
           l_gga= xcpot%is_gga()
           l_test = .TRUE.                  ! only checking, dont use new parameters
-          CALL chkmt(atoms,input,vacuum,cell,oneD,l_gga,noel,l_test, kmax1,dtild,dvac1,lmax1,jri1,rmt1,dx1)
+          CALL chkmt(atoms,input,vacuum,cell,oneD,l_test,l_gga,noel, kmax1,dtild,dvac1,lmax1,jri1,rmt1,dx1)
 
           WRITE (6,FMT=8180) cell%volint
 8180      FORMAT (13x,' volume of interstitial region=',f12.6)

init/postprocessInput.F90

@@ -38,6 +38,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
   USE m_nocoInputCheck
   USE m_kpoints   
   USE m_types_forcetheo_extended
+  USE m_relaxio
 
   IMPLICIT NONE
 
@@ -319,15 +320,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
      END DO
 
 
-     ! Check muffin tin radii
 
-     ALLOCATE (jri1(atoms%ntype), lmax1(atoms%ntype))
-     ALLOCATE (rmt1(atoms%ntype), dx1(atoms%ntype))
-     l_test = .TRUE. ! only checking, dont use new parameters
-     l_gga=xcpot%is_gga()
-     CALL chkmt(atoms,input,vacuum,cell,oneD,l_gga,noel,l_test,&
-                kmax1,dtild1,dvac1,lmax1,jri1,rmt1,dx1)
-     DEALLOCATE (jri1,lmax1,rmt1,dx1)
 
      ! Dimensioning of lattice harmonics
 
@@ -468,6 +461,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
      IF (atoms%n_u.GT.0) THEN
         CALL d_wigner(sym%nop,sym%mrot,cell%bmat,3,sym%d_wgn)
      END IF
+
      IF (.NOT.oneD%odd%d1) THEN
         CALL mapatom(sym,atoms,cell,input,noco)
         oneD%ngopr1(1:atoms%nat) = atoms%ngopr(1:atoms%nat)
@@ -477,6 +471,13 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
         CALL mapatom(sym,atoms,cell,input,noco)
         CALL od_mapatom(oneD,atoms,sym,cell)
      END IF
+    ! Check muffin tin radii
+
+     l_test = .TRUE. ! only checking, dont use new parameters
+     CALL chkmt(atoms,input,vacuum,cell,oneD,l_test)
+
+     !adjust positions by displacements
+     CALL apply_displacements(cell,input,vacuum,oneD,sym,noco,atoms)
 
      !Calculate kpoint in the full BZ
      IF (kpts%l_gamma.and. banddos%ndir .eq. 0.and.kpts%specificationType==2) THEN
@@ -551,7 +552,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
      END IF !(mpi%irank.EQ.0)
   END IF
 
-  ! 
+ 
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!! End of input postprocessing (calculate missing parameters)

inpgen/set_inp.f90

@@ -201,8 +201,8 @@
       a1(:) = cell%amat(:,1) ; a2(:) = cell%amat(:,2) ; a3(:) = cell%amat(:,3) 
 
       CALL chkmt(&
-     &           atoms,input,vacuum,cell,oneD,&
-     &           l_gga,noel,l_test,&
+     &           atoms,input,vacuum,cell,oneD,l_test,&
+     &           l_gga,noel,&
      &           kmax,dtild,vacuum%dvac,atoms%lmax,atoms%jri,atoms%rmt,atoms%dx)
 
 ! --> read in (possibly) atomic info
@@ -236,8 +236,8 @@
       rmtTemp = 999.0
       l_test = .true.
       CALL chkmt(&
-     &           atoms,input,vacuum,cell,oneD,&
-     &           l_gga,noel,l_test,&
+     &           atoms,input,vacuum,cell,oneD,l_test,&
+     &           l_gga,noel,&
      &           kmax0,dtild0,dvac0,lmax0,jri0,rmtTemp,dx0)
 
       IF ( ANY(atoms%nlo(:).NE.0) ) THEN

io/r_inpXML.F90

@@ -704,6 +704,10 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
      
       xPathA = '/fleurInput/calculationSetup/fields'
       numberNodes = xmlGetNumberOfNodes(xPathA)
+      field%b_field=0.0
+      field%l_b_field=.FALSE.
+      field%efield%sigma=0.0
+      
 
       IF (numberNodes.EQ.1) THEN
          IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@b_field')>0) THEN
@@ -725,7 +729,8 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
             WRITE(xPathB,"(a,a,i0,a)") TRIM(ADJUSTL(xpathA)),'/shape[',i,']'
             field%efield%shapes(i)=TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
          ENDDO
-         
+      ELSE
+         ALLOCATE(field%efield%shapes(0))
       END IF
 
       ! Read in optional energy parameter limits
@@ -2139,12 +2144,13 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
          CALL inpnoco(atoms,input,vacuum,noco)
       END IF
 
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!! End of non-XML input
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       !CALL xmlFreeResources()
-
+      
       !WRITE(*,*) 'Reading of inp.xml file finished'
 
       DEALLOCATE(speciesNLO)

io/relax_io.F90

@@ -8,7 +8,7 @@ MODULE m_relaxio
   USE m_judft
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: read_relax,write_relax
+  PUBLIC :: read_relax,write_relax,apply_displacements,read_displacements
 CONTAINS
   SUBROUTINE write_relax(positions,forces,energies,displace)
     REAL,INTENT(in):: positions(:,:,:)
@@ -19,16 +19,17 @@ CONTAINS
     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")
+    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
+    OPEN(765,file="relax.xml",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>'
+       WRITE(765,"(a,3(f15.10,1x),a)") &
+            '    <displace>',displace(:,n),'</displace>'
     END DO
     WRITE(765,"(a)") '  </displacements>'
     
@@ -37,10 +38,8 @@ CONTAINS
     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>'
+          WRITE(765,"(a,6(f15.10,1x),a)") &
+               '      <posforce>',positions(:,n,step),forces(:,n,step),'</posforce>'
        END DO
        WRITE(765,"(a)") '    </step>'
     ENDDO
@@ -52,61 +51,144 @@ CONTAINS
   SUBROUTINE read_relax(positions,forces,energies)
     USE m_xmlIntWrapFort 
     USE m_calculator
-    REAL,INTENT(OUT),ALLOCATABLE:: positions(:,:,:)