Feature/print stress t25

cmake/CMakeLists.txt

@@ -144,16 +144,16 @@ set(_FILES
     params paths pbstuf pdb phipsi picalc piorbs pistuf pitors pme
     pmestuf pmpb polar polgrp polopt polpcg polpot poltcg polymer potent
     potfit predict pressure prmkey promo prtarc prtdyn prterr prtfrc
-    prtint prtmol2 prtpdb prtprm prtseq prtuind prtvel prtxyz ptable
+    prtint prtmol2 prtpdb prtprm prtseq prtstres prtuind prtvel prtxyz ptable
     qmstuf qrsolve quatfit random rattle readcart readdcd readdyn
     readgau readgdma readint readmbis readmol readmol2 readpdb readprm
     readseq readxyz refer repel replica reppot resdue respa restrn rgddyn
     rgdstep richmond rigid ring rings rmsfit rotbnd rotlist rotpole
     rxnfld rxnpot scales sdstep search sequen server setprm shakeup
     shapes shunt sigmoid simplex sizes sktstuf socket solpot solute sort
-    square stodyn strbnd strtor suffix surface surfatom switch syntrn
-    tarray tcgstuf temper tettor tettors titles tncg torphase torpot
-    torque tors torsions tortor tree trimtext tritor tritors unionball
+    square stodyn strbnd stresinit strtor strvar suffix surface surfatom
+    switch syntrn tarray tcgstuf temper tettor tettors titles tncg torphase
+    torpot torque tors torsions tortor tree trimtext tritor tritors unionball
     unitcell units uprior urey urypot usage valfit vdw vdwpot verlet
     version vibs virial volume warp xtals xyzatm zatom zclose zcoord
 )

source/beeman.f

@@ -159,6 +159,7 @@ subroutine beeman (istep,dt)
 c
 c     compute statistics and save trajectory for this step
 c
+      call prtstres (istep,stress)
       call mdstat (istep,dt,etot,epot,eksum,temp,pres)
       call mdsave (istep,dt,epot,eksum)
       call mdrest (istep)

source/dynamic.f

@@ -227,6 +227,10 @@ program dynamic
 c
       call mdinit (dt)
 c
+c     check for keywords that alter stress tensor output
+c
+      call stresinit (dt)
+c
 c     print out a header line for the dynamics computation
 c
       if (integrate .eq. 'VERLET') then

source/prtstres.f

@@ -0,0 +1,72 @@
+c
+c
+c     ###################################################
+c     ##                                               ##
+c     ##                                               ##
+c     ###################################################
+c
+c     ################################################################
+c     ##                                                            ##
+c     ##  subroutine prtstres  --  output of stress tensor          ##
+c     ##                                                            ##
+c     ################################################################
+c
+c
+c     "prtstres" writes stress tensor components to an output file
+c
+c
+      subroutine prtstres (istep,stress)
+      use bound
+      use boxes
+      use files
+      use inform
+      use iounit
+      use output
+      use strvar
+      use titles
+      implicit none
+      integer istep
+      integer istr
+      integer modstre
+      integer freeunit
+      real*8 stress(3,3)
+      logical exist
+      character*240 strfile
+c
+c
+c     only necessary if user requests stress tensor output
+c
+      if (.not. stresav)  return
+c
+c     check number of steps between trajectory file saves
+c
+      modstre = mod(istep,istress)
+      if (modstre .ne. 0)  return
+c
+c     open and format the output file 
+c
+      istr = freeunit ()
+      strfile = filename(1:leng)
+      call suffix (strfile,'str','old')
+      inquire (file=strfile,exist=exist)
+      if (exist) then
+         call openend (istr,strfile)
+      else
+         open (unit=istr,file=strfile,status='new')
+         write (istr, 10)
+   10    format ('MD STEP',4x,
+     &           'STRESS(1,1)',4x,'STRESS(2,1)',4x,'STRESS(3,1)',4x,
+     &           'STRESS(1,2)',4x,'STRESS(2,2)',4x,'STRESS(3,2)',4x,
+     &           'STRESS(1,3)',4x,'STRESS(2,3)',4x,'STRESS(3,3)')
+      end if
+c
+c     write the value of the stress tensor components
+c
+      write (istr, 20) istep,stress
+   20 format (i7,9f15.4)
+c
+c     close the stress tensor output file
+c      
+      close (unit=istr)
+      return
+      end

source/stresinit.f

@@ -0,0 +1,52 @@
+c
+c
+c     ###################################################
+c     ##                                               ##
+c     ##                                               ##
+c     ###################################################
+c
+c     ################################################################
+c     ##                                                            ##
+c     ##  subroutine stresinit  --  set stress tenor output         ##
+c     ##                                                            ##
+c     ################################################################
+c
+c
+c     "stresinit" checks if user requests outputting the stress tensor  
+c     and if so, how often.
+c
+c
+c
+      subroutine stresinit (dt)
+      use keys
+      use strvar
+      implicit none
+      integer i,next
+      real*8 dt
+      character*20 keyword
+      character*240 record
+      character*240 string
+c
+c
+c     set default values for stress tensor output variables
+c
+      stresav = .false.
+      istress = 10
+c
+c     search keywords for stress output parameters
+c
+      do i = 1, nkey
+         next = 1
+         record = keyline(i)
+         call gettext (record,keyword,next)
+         call upcase (keyword)
+         string = record(next:240)
+         if (keyword(1:12) .eq. 'SAVE-STRESS ') then
+            stresav = .true.
+         else if (keyword(1:12) .eq. 'STRESS-FREQ ') then
+            read (string,*,err=10,end=10) istress
+         end if
+   10    continue
+      end do
+      return
+      end

source/strvar.f

@@ -0,0 +1,25 @@
+c
+c
+c     ###################################################
+c     ##                                               ##
+c     ##                                               ##
+c     ###################################################
+c
+c     #################################################################
+c     ##                                                             ##
+c     ##  module strvar -- stress tensor output control parameters   ##
+c     ##                                                             ##
+c     #################################################################
+c
+c
+c     stresav   logical flag to save stress tensor components
+c     istress   steps between stress tensor component prints   
+c
+c
+      module strvar
+      implicit none
+      integer istress
+      logical stresav
+      save
+      end
+

source/verlet.f

@@ -23,15 +23,17 @@ subroutine verlet (istep,dt)
       use ielscf
       use moldyn
       use polar
+      use strvar
       use units
       use usage
       implicit none
       integer i,j,k
       integer istep
       real*8 dt,dt_2
       real*8 etot,epot
-      real*8 eksum,term
+      real*8 eksum
       real*8 temp,pres
+      real*8 term
       real*8 ekin(3,3)
       real*8 stress(3,3)
       real*8, allocatable :: xold(:)
@@ -142,6 +144,7 @@ subroutine verlet (istep,dt)
 c
 c     compute statistics and save trajectory for this step
 c
+      call prtstres (istep,stress)
       call mdstat (istep,dt,etot,epot,eksum,temp,pres)
       call mdsave (istep,dt,epot,eksum)
       call mdrest (istep)