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