1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  subroutine mdsave  --  save trajectory and restart files  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "mdsave" writes molecular dynamics trajectory snapshots and
16c     auxiliary files with velocity, force or induced dipole data;
17c     also checks for user requested termination of a simulation
18c
19c
20      subroutine mdsave (istep,dt,epot,eksum)
21      use atomid
22      use atoms
23      use bound
24      use boxes
25      use deriv
26      use files
27      use group
28      use inform
29      use iounit
30      use mdstuf
31      use moldyn
32      use mpole
33      use output
34      use polar
35      use potent
36      use rgddyn
37      use socket
38      use titles
39      use units
40      implicit none
41      integer i,j,k,istep
42      integer ixyz,iind
43      integer ivel,ifrc
44      integer iend,isave,lext
45      integer freeunit,trimtext
46      integer modsave
47      real*8 dt,pico
48      real*8 epot,eksum
49      logical exist
50      character*7 ext
51      character*240 endfile
52      character*240 xyzfile
53      character*240 velfile
54      character*240 frcfile
55      character*240 indfile
56c
57c
58c     send data via external socket communication if desired
59c
60      if (.not.sktstart .or. use_socket)  call sktdyn (istep,dt,epot)
61c
62c     check number of steps between trajectory file saves
63c
64      modsave = mod(istep,iwrite)
65      if (modsave .ne. 0)  return
66c
67c     get the sequence number of the current trajectory frame
68c
69      isave = nprior + istep/iwrite
70      lext = 3
71      call numeral (isave,ext,lext)
72c
73c     print header for the instantaneous values at current step
74c
75      pico = dble(istep) * dt
76      write (iout,10)  istep
77   10 format (/,' Instantaneous Values for Frame Saved at',
78     &           i10,' Dynamics Steps')
79c
80c     print the current time, potential and kinetic energies
81c
82      if (digits .ge. 8) then
83         write (iout,20)  pico
84   20    format (/,' Current Time',8x,f19.8,' Picosecond')
85         write (iout,30)  epot
86   30    format (' Current Potential',3x,f19.8,' Kcal/mole')
87         write (iout,40)  eksum
88   40    format (' Current Kinetic',5x,f19.8,' Kcal/mole')
89      else if (digits .ge. 6) then
90         write (iout,50)  pico
91   50    format (/,' Current Time',8x,f17.6,' Picosecond')
92         write (iout,60)  epot
93   60    format (' Current Potential',3x,f17.6,' Kcal/mole')
94         write (iout,70)  eksum
95   70    format (' Current Kinetic',5x,f17.6,' Kcal/mole')
96      else
97         write (iout,80)  pico
98   80    format (/,' Current Time',8x,f15.4,' Picosecond')
99         write (iout,90)  epot
100   90    format (' Current Potential',3x,f15.4,' Kcal/mole')
101         write (iout,100)  eksum
102  100    format (' Current Kinetic',5x,f15.4,' Kcal/mole')
103      end if
104c
105c     print the values of the lattice lengths and angles
106c
107      if (use_bounds) then
108         if (digits .le. 6) then
109            write (iout,110)  xbox,ybox,zbox
110  110       format (' Lattice Lengths',6x,3f14.6)
111            write (iout,120)  alpha,beta,gamma
112  120       format (' Lattice Angles',7x,3f14.6)
113         else if (digits .le. 8) then
114            write (iout,130)  xbox,ybox,zbox
115  130       format (' Lattice Lengths',6x,3f16.8)
116            write (iout,140)  alpha,beta,gamma
117  140       format (' Lattice Angles',7x,3f16.8)
118         else
119            write (iout,150)  xbox,ybox,zbox
120  150       format (' Lattice Lengths',6x,3f18.10)
121            write (iout,160)  alpha,beta,gamma
122  160       format (' Lattice Angles',7x,3f18.10)
123         end if
124      end if
125c
126c     save coordinates to an archive or numbered structure file
127c
128      ixyz = freeunit ()
129      if (archive) then
130         xyzfile = filename(1:leng)
131         call suffix (xyzfile,'arc','old')
132         inquire (file=xyzfile,exist=exist)
133         if (exist) then
134            call openend (ixyz,xyzfile)
135         else
136            open (unit=ixyz,file=xyzfile,status='new')
137         end if
138      else
139         xyzfile = filename(1:leng)//'.'//ext(1:lext)
140         call version (xyzfile,'new')
141         open (unit=ixyz,file=xyzfile,status='new')
142      end if
143      if (use_bounds)  call bounds
144      call prtxyz (ixyz)
145      close (unit=ixyz)
146      write (iout,170)  isave
147  170 format (' Frame Number',13x,i10)
148      write (iout,180)  xyzfile(1:trimtext(xyzfile))
149  180 format (' Coordinate File',13x,a)
150c
151c     update the information needed to restart the trajectory
152c
153      call prtdyn
154c
155c     save the velocity vector components at the current step
156c
157      if (velsave) then
158         ivel = freeunit ()
159         if (archive) then
160            velfile = filename(1:leng)
161            call suffix (velfile,'vel','old')
162            inquire (file=velfile,exist=exist)
163            if (exist) then
164               call openend (ivel,velfile)
165            else
166               open (unit=ivel,file=velfile,status='new')
167            end if
168         else
169            velfile = filename(1:leng)//'.'//ext(1:lext)//'v'
170            call version (velfile,'new')
171            open (unit=ivel,file=velfile,status='new')
172         end if
173         if (integrate .eq. 'RIGIDBODY') then
174            write (ivel,190)  ngrp,title(1:ltitle)
175  190       format (i6,2x,a)
176            do i = 1, ngrp
177               write (ivel,200)  i,(vcm(j,i),j=1,3)
178  200          format (i6,3x,d13.6,3x,d13.6,3x,d13.6)
179               write (ivel,210)  i,(wcm(j,i),j=1,3)
180  210          format (i6,3x,d13.6,3x,d13.6,3x,d13.6)
181            end do
182         else
183            write (ivel,220)  n,title(1:ltitle)
184  220       format (i6,2x,a)
185            do i = 1, n
186               write (ivel,230)  i,name(i),(v(j,i),j=1,3)
187  230          format (i6,2x,a3,3x,d13.6,3x,d13.6,3x,d13.6)
188            end do
189         end if
190         close (unit=ivel)
191         write (iout,240)  velfile(1:trimtext(velfile))
192  240    format (' Velocity File',15x,a)
193      end if
194c
195c     save the force vector components for the current step,
196c     only correct for single time step Cartesian integrators
197c
198      if (frcsave .and. integrate.ne.'RIGIDBODY'
199     &       .and. integrate.ne.'RESPA') then
200         ifrc = freeunit ()
201         if (archive) then
202            frcfile = filename(1:leng)
203            call suffix (frcfile,'frc','old')
204            inquire (file=frcfile,exist=exist)
205            if (exist) then
206               call openend (ifrc,frcfile)
207            else
208               open (unit=ifrc,file=frcfile,status='new')
209            end if
210         else
211            frcfile = filename(1:leng)//'.'//ext(1:lext)//'f'
212            call version (frcfile,'new')
213            open (unit=ifrc,file=frcfile,status='new')
214         end if
215         write (ifrc,250)  n,title(1:ltitle)
216  250    format (i6,2x,a)
217         do i = 1, n
218            write (ifrc,260)  i,name(i),(-desum(j,i),j=1,3)
219  260       format (i6,2x,a3,3x,d13.6,3x,d13.6,3x,d13.6)
220         end do
221         close (unit=ifrc)
222         write (iout,270)  frcfile(1:trimtext(frcfile))
223  270    format (' Force Vector File',11x,a)
224      end if
225c
226c     save the current induced dipole moment at each site
227c
228      if (uindsave .and. use_polar) then
229         iind = freeunit ()
230         if (archive) then
231            indfile = filename(1:leng)
232            call suffix (indfile,'uind','old')
233            inquire (file=indfile,exist=exist)
234            if (exist) then
235               call openend (iind,indfile)
236            else
237               open (unit=iind,file=indfile,status='new')
238            end if
239         else
240            indfile = filename(1:leng)//'.'//ext(1:lext)//'u'
241            call version (indfile,'new')
242            open (unit=iind,file=indfile,status='new')
243         end if
244         write (iind,280)  n,title(1:ltitle)
245  280    format (i6,2x,a)
246         do i = 1, npole
247            if (polarity(i) .ne. 0.0d0) then
248               k = ipole(i)
249               write (iind,290)  k,name(k),(debye*uind(j,i),j=1,3)
250  290          format (i6,2x,a3,3f12.6)
251            end if
252         end do
253         close (unit=iind)
254         write (iout,300)  indfile(1:trimtext(indfile))
255  300    format (' Induced Dipole File',9x,a)
256      end if
257c
258c     test for requested termination of the dynamics calculation
259c
260      endfile = 'tinker.end'
261      inquire (file=endfile,exist=exist)
262      if (.not. exist) then
263         endfile = filename(1:leng)//'.end'
264         inquire (file=endfile,exist=exist)
265         if (exist) then
266            iend = freeunit ()
267            open (unit=iend,file=endfile,status='old')
268            close (unit=iend,status='delete')
269         end if
270      end if
271      if (exist) then
272         write (iout,310)
273  310    format (/,' MDSAVE  --  Dynamics Calculation Ending',
274     &              ' due to User Request')
275         call fatal
276      end if
277c
278c     skip an extra line to keep the output formating neat
279c
280      modsave = mod(istep,iprint)
281      if (verbose .and. modsave.ne.0) then
282         write (iout,320)
283  320    format ()
284      end if
285      return
286      end
287