From 7cbf956031f8f665c329a301d608fed88b536e1c Mon Sep 17 00:00:00 2001 From: Daniel Relix Date: Wed, 4 Jun 2025 17:19:29 -0500 Subject: [PATCH 1/7] Add stress tensor output module and subroutines --- source/prtstres.f | 71 ++++++++++++++++++++++++++++++++++++++++++++++ source/stresinit.f | 59 ++++++++++++++++++++++++++++++++++++++ source/strvar.f | 27 ++++++++++++++++++ 3 files changed, 157 insertions(+) create mode 100644 source/prtstres.f create mode 100644 source/stresinit.f create mode 100644 source/strvar.f diff --git a/source/prtstres.f b/source/prtstres.f new file mode 100644 index 00000000..60cad3dd --- /dev/null +++ b/source/prtstres.f @@ -0,0 +1,71 @@ +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 + 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 \ No newline at end of file diff --git a/source/stresinit.f b/source/stresinit.f new file mode 100644 index 00000000..d58048f2 --- /dev/null +++ b/source/stresinit.f @@ -0,0 +1,59 @@ +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 NOTE: This is the minimum amount of code required for functionality. +c A future "stresinit" update could add logic to verify that correct +c ensemble (NVT) and thermostat (Velocity-Scaling) are selected. +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. + stresfrq = 10.0d0 +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) stresfrq + end if + 10 continue + end do + if (stresav) then + stresfrq = 0.001d0 * stresfrq + istress = nint(stresfrq/dt) + end if + return + end diff --git a/source/strvar.f b/source/strvar.f new file mode 100644 index 00000000..5ba80ef7 --- /dev/null +++ b/source/strvar.f @@ -0,0 +1,27 @@ +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 stresfrq time (in femtoseconds) between stress tensor prints +c istress steps between stress tensor component prints +c +c + module strvar + implicit none + integer istress + real*8 stresfrq + logical stresav + save + end + From 15abcfbb19c484217edb7552d5cec4cf3e0c12fd Mon Sep 17 00:00:00 2001 From: Daniel Relix Date: Wed, 4 Jun 2025 17:24:57 -0500 Subject: [PATCH 2/7] Modify dynamic.f and verlet.f to call stress output subroutines --- source/dynamic.f | 4 ++++ source/verlet.f | 5 ++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/source/dynamic.f b/source/dynamic.f index 491bbc67..c8e80c3d 100644 --- a/source/dynamic.f +++ b/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 diff --git a/source/verlet.f b/source/verlet.f index 51de0f16..cbb81c08 100644 --- a/source/verlet.f +++ b/source/verlet.f @@ -23,6 +23,7 @@ subroutine verlet (istep,dt) use ielscf use moldyn use polar + use strvar use units use usage implicit none @@ -30,8 +31,9 @@ subroutine verlet (istep,dt) 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) From c9491a284ebb8f8aa79933b4c5088e10930a3874 Mon Sep 17 00:00:00 2001 From: Daniel Relix Date: Wed, 4 Jun 2025 21:04:59 -0500 Subject: [PATCH 3/7] Update CMakeLists.txt to include new source files --- cmake/CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index e65487cb..8deb9c1c 100644 --- a/cmake/CMakeLists.txt +++ b/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 ) From 8dbe3faaeea893b3bc135e059c6aa332ad1f4d9d Mon Sep 17 00:00:00 2001 From: Daniel Relix Date: Sun, 28 Sep 2025 18:44:33 -0500 Subject: [PATCH 4/7] Update beeman.f to support stress tensor output --- source/beeman.f | 1 + 1 file changed, 1 insertion(+) diff --git a/source/beeman.f b/source/beeman.f index af22bf87..79fe49d8 100644 --- a/source/beeman.f +++ b/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) From f4c96fa97c30f8c77abfae96a6a3ea683469e6df Mon Sep 17 00:00:00 2001 From: Daniel Relix Date: Sun, 28 Sep 2025 19:04:14 -0500 Subject: [PATCH 5/7] Update prtstres.f for stress output handling --- source/prtstres.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/prtstres.f b/source/prtstres.f index 60cad3dd..1db0faba 100644 --- a/source/prtstres.f +++ b/source/prtstres.f @@ -28,6 +28,7 @@ subroutine prtstres (istep,stress) integer istep integer istr integer modstre + integer freeunit real*8 stress(3,3) logical exist character*240 strfile @@ -68,4 +69,4 @@ subroutine prtstres (istep,stress) c close (unit=istr) return - end \ No newline at end of file + end From 58b315cd3f93f9f25ceb6132dd07cd7fe54a23d0 Mon Sep 17 00:00:00 2001 From: Daniel Relix Date: Wed, 12 Nov 2025 17:24:22 -0600 Subject: [PATCH 6/7] Remove uneeded variables --- source/stresinit.f | 11 ++--------- source/strvar.f | 1 - 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/source/stresinit.f b/source/stresinit.f index d58048f2..8fdcfd58 100644 --- a/source/stresinit.f +++ b/source/stresinit.f @@ -15,9 +15,6 @@ c "stresinit" checks if user requests outputting the stress tensor c and if so, how often. c -c NOTE: This is the minimum amount of code required for functionality. -c A future "stresinit" update could add logic to verify that correct -c ensemble (NVT) and thermostat (Velocity-Scaling) are selected. c c subroutine stresinit (dt) @@ -34,7 +31,7 @@ subroutine stresinit (dt) c set default values for stress tensor output variables c stresav = .false. - stresfrq = 10.0d0 + istress = 10 c c search keywords for stress output parameters c @@ -47,13 +44,9 @@ subroutine stresinit (dt) 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) stresfrq + read (string,*,err=10,end=10) istress end if 10 continue end do - if (stresav) then - stresfrq = 0.001d0 * stresfrq - istress = nint(stresfrq/dt) - end if return end diff --git a/source/strvar.f b/source/strvar.f index 5ba80ef7..cf03311d 100644 --- a/source/strvar.f +++ b/source/strvar.f @@ -20,7 +20,6 @@ module strvar implicit none integer istress - real*8 stresfrq logical stresav save end From 788bdf0a8de254d8f81f10e006b556ee94ff6fad Mon Sep 17 00:00:00 2001 From: Daniel R <115332319+d-rel82@users.noreply.github.com> Date: Thu, 13 Nov 2025 13:16:51 -0600 Subject: [PATCH 7/7] Update strvar.f --- source/strvar.f | 1 - 1 file changed, 1 deletion(-) diff --git a/source/strvar.f b/source/strvar.f index cf03311d..b31c4e4b 100644 --- a/source/strvar.f +++ b/source/strvar.f @@ -13,7 +13,6 @@ c c c stresav logical flag to save stress tensor components -c stresfrq time (in femtoseconds) between stress tensor prints c istress steps between stress tensor component prints c c