diff --git a/GPU/source/MOD_mdstuf.f b/GPU/source/MOD_mdstuf.f
index de36c1a8..5a3edefc 100644
--- a/GPU/source/MOD_mdstuf.f
+++ b/GPU/source/MOD_mdstuf.f
@@ -16,6 +16,7 @@
 c     dorest      logical flag to remove center of mass inertia
 c     velsave     logical flag to save velocity vector components
 c     frcsave     logical flag to save force vector components
+c     stresave    logical flag to save stress tensor components
 c     uindsave    logical flag to save induced atomic dipoles
 c     integrate   type of molecular dynamics integration algorithm
 c
@@ -28,6 +29,7 @@ module mdstuf
       logical dorest
       logical velsave
       logical frcsave
+      logical stresave
       logical uindsave
       logical mts
       character*20 integrate
diff --git a/GPU/source/mdinit.f b/GPU/source/mdinit.f
index 2cc626df..10b6a27a 100644
--- a/GPU/source/mdinit.f
+++ b/GPU/source/mdinit.f
@@ -212,6 +212,13 @@ function maxwell (mass,temper)
           frcsave = .true.
         case ('SAVE-INDUCED')
           uindsave = .true.
+        case ('SAVE-STRESS')
+          stresave = .true.
+          use_virial=.true.
+          verbose = .true.
+          write(*,*)
+     &    'SAVE-STRESS: The units of the stress tensor '// 
+     &    'components are Atmospheres'
         case ('THERMOSTAT')
           call getword (record,thermostat,next)
           call upcase (thermostat)
diff --git a/GPU/source/mdstat.f b/GPU/source/mdstat.f
index 69d9182c..c7a71245 100644
--- a/GPU/source/mdstat.f
+++ b/GPU/source/mdstat.f
@@ -227,6 +227,15 @@ subroutine mdstat (istep,dt,etot,epot,ekin,temp,pres)
                    write(iout,'(A)',advance="no") 
      &                         "       Pres"
                  endif
+                 if(stresave) then
+                   write(iout,'(A)',advance="no")
+     &                         "     Stress(xx)"
+     &                       //"     Stress(yy)"
+     &                       //"     Stress(zz)"
+     &                       //"     Stress(xy)"
+     &                       //"     Stress(yz)"
+     &                       //"     Stress(xz)"  
+                 endif      
 c                 if(isobaric) then
 c                   write(iout,'(A)',advance="no") 
 c     &                         "     Density"
@@ -243,6 +252,11 @@ subroutine mdstat (istep,dt,etot,epot,ekin,temp,pres)
                if(use_virial) then
                  write(iout,'(f11.2)',advance="no") pres
                endif
+               if(stresave) then
+                 write(iout,'(6f15.4)',advance="no") 
+     &                 stress(1,1),stress(2,2),stress(3,3),
+     &                 stress(1,2),stress(2,3),stress(1,3)
+               endif
                !if(isobaric) then
                !  write(iout,'(f12.4,f12.2)',advance="no") dens,volbox
                !endif
diff --git a/GPU/source/pressure.f b/GPU/source/pressure.f
index 47cb2176..52b026bd 100644
--- a/GPU/source/pressure.f
+++ b/GPU/source/pressure.f
@@ -21,6 +21,7 @@ subroutine pressure (dt,ekin,pres,stress,istep)
       use bound
       use boxes
       use domdec
+      use mdstuf    ,only:stresave
       use tinheader ,only:ti_p,re_p
       use units
       use virial
@@ -34,14 +35,13 @@ subroutine pressure (dt,ekin,pres,stress,istep)
       real(r_p) stres1,stres2,stres3
 c
 c     only necessary if periodic boundaries are in use
-c     and isobaric simulation
 c
       if (.not.(use_bounds.and.use_virial))  return
 c
 c     calculate the stress tensor for anisotropic systems
 c
       factor = prescon / volbox
-      if (anisotrop) then
+      if (anisotrop.or.stresave) then
 !$acc parallel loop collapse(2) default(present) async
         do i = 1, 3
            do j = 1, 3