diff --git a/force/geo.f90 b/force/geo.f90
index 086c072961d5d6a81330b9a84c3f4e405080ccc5..be66d194076bd9b6e25409017f453542d594c1cf 100644
--- a/force/geo.f90
+++ b/force/geo.f90
@@ -42,6 +42,7 @@ CONTAINS
     USE m_rinpXML
     USE m_winpXML
     USE m_init_wannier_defaults
+    USE m_xsf_io
     IMPLICIT NONE
     TYPE(t_oneD),INTENT(IN)   :: oneD
     TYPE(t_cell),INTENT(IN)   :: cell
@@ -101,6 +102,7 @@ CONTAINS
     CHARACTER(LEN=20)             :: filename
     REAL                          :: a1_temp(3),a2_temp(3),a3_temp(3)
     REAL                          :: scale_temp, dtild_temp
+    REAL                          :: forceAllAtoms(3,atoms%nat)
     CLASS(t_forcetheo),ALLOCATABLE:: forcetheo
     input=input_in
     atoms_new=atoms
@@ -135,9 +137,33 @@ CONTAINS
     ENDDO
 
     istep = 1
-    CALL bfgs(atoms%ntype,istep,istep0,forcetot,&
-         &          zat,input%xa,input%thetad,input%epsdisp,input%epsforce,tote,&
-         &          xold,y,h,tau0, lconv)
+    CALL bfgs(atoms%ntype,istep,istep0,forcetot,zat,input%xa,input%thetad,input%epsdisp,&
+              input%epsforce,tote,xold,y,h,tau0,lconv)
+
+    !write out struct_force.xsf file
+    forceAllAtoms = 0.0
+    na = 0
+    DO itype=1,atoms%ntype
+       forcetot(:,itype)=MATMUL(cell%bmat,forcetot(:,itype))/tpi_const ! to inner coordinates
+       DO ieq = 1,atoms%neq(itype)
+          na = na + 1
+          jop = sym%invtab(atoms%ngopr(na))
+          IF (oneD%odi%d1) jop = oneD%ods%ngopr(na)
+          DO i = 1,3
+             DO j = 1,3
+                IF (.NOT.oneD%odi%d1) THEN
+                   forceAllAtoms(i,na) = forceAllAtoms(i,na) + sym%mrot(i,j,jop) * forcetot(j,itype)
+                ELSE
+                   forceAllAtoms(i,na) = forceAllAtoms(i,na) + oneD%ods%mrot(i,j,jop) * forcetot(j,itype)
+                END IF
+             END DO
+          END DO
+          forceAllAtoms(:,na) = MATMUL(cell%amat,forceAllAtoms(:,na)) ! to external coordinates
+       END DO
+    END DO
+    OPEN (55,file="struct_force.xsf",status='replace')
+    CALL xsf_WRITE_atoms(55,atoms,input%film,.false.,cell%amat,forceAllAtoms)
+    CLOSE (55)
 
     IF (lconv) THEN
        WRITE (6,'(a)') "Des woars!"
diff --git a/io/xsf_io.f90 b/io/xsf_io.f90
index c8b8be93f74f319ca9640293eeb41ea3f70df4bd..cb1600840c768e9ef4c017a757dcd20b9cb3a612 100644
--- a/io/xsf_io.f90
+++ b/io/xsf_io.f90
@@ -15,7 +15,7 @@ MODULE m_xsf_io
 
 CONTAINS
   !<-- S:S: xsf_WRITE_atoms(fileno,film,amat,neq(:ntype),zatom(:ntype),pos)
-  SUBROUTINE xsf_WRITE_atoms(fileno,atoms,film,od,amat)
+  SUBROUTINE xsf_WRITE_atoms(fileno,atoms,film,od,amat,forceAllAtoms)
     !-----------------------------------------------
     !     Writes the crystal dimensions&atomic positions
     !           (last modified: 2004-00-00) D. Wortmann
@@ -27,6 +27,7 @@ CONTAINS
     LOGICAL,INTENT(IN)      :: film
     LOGICAL,INTENT(IN)     :: od
     REAL,INTENT(IN)        :: amat(3,3)
+    REAL, OPTIONAL, INTENT(IN) :: forceAllAtoms(3,atoms%nat)
     !>
     !<-- Locals
     INTEGER             :: n,nn,na
@@ -52,8 +53,11 @@ CONTAINS
     na = 1
     DO n = 1,SIZE(atoms%neq)
        DO nn = 1,atoms%neq(n)
-          WRITE(fileno,'(i4,2x,3(f0.7,1x))') NINT(atoms%zatom(n)),&
-               &            atoms%pos(:,na)*a0
+          IF (PRESENT(forceAllAtoms)) THEN
+             WRITE(fileno,'(i4,2x,6(f0.7,1x))') NINT(atoms%zatom(n)),atoms%pos(:,na)*a0,forceAllAtoms(:,na)/a0
+          ELSE
+             WRITE(fileno,'(i4,2x,3(f0.7,1x))') NINT(atoms%zatom(n)),atoms%pos(:,na)*a0
+          END IF
           na=na+1
        ENDDO
     ENDDO