program ham implicit none integer*4 NW,i,NEQUI,NENS,NPROD,i1,i2,i3,i4, 1istart,iend,j,ist,ien real*8 xmin,xmax,TEMP,NFS,NPS,PRES,x,dx,b,MTOL character*10 so,sop character*80 s80,dynamic,minimize logical lex,LPER,lpc,ltinker5 c dynamic='dynamic' minimize='minimize' ltinker5=.false. MTOL=1.0d0 NEQUI=1000 NPROD=10000 NFS=1.0d0 NENS=2 TEMP=298.0d0 PRES=1.0d0 XMIN=0.0d0 XMAX=0.0d0 NW=50 LPER=.false. lpc=.false. write(6,*) write(6,600) 600 format(' Makes the driving script for WHAM analysis in Tinker',/, 1 ' ',/, 1 ' Petr Bour, Academy of Sciences, 2010 ',/, 1 ' ',/, 1 ' Input: HAM.OPT',/, 1 ' KEY.TXT (the .key file torso)',/, 1 ' ham.xyz.or (for MD)',/, 1 ' ham.prm (for MD)',/, 1 ' ',/, 1 ' Output: goham',/,/) open(9,file='HAM.OPT') 1 read(9,900,end=99,err=99)so 900 format(a10) if(so(1:4).eq.'DYNA')read(9,'(a)')dynamic if(so(1:4).eq.'LPER')read(9,*)LPER if(so(1:4).eq.'TINK')read(9,*)ltinker5 if(so(1:4).eq.'MINI')read(9,'(a)')minimize if(so(1:4).eq.'MTOL')read(9,*)MTOL if(so(1:5).eq.'NEQUI')read(9,*)NEQUI if(so(1:4).eq.'NENS')read(9,*)NENS if(so(1:3).eq.'NFS')read(9,*)NFS if(so(1:4).eq.'NINT')read(9,*) if(so(1:5).eq.'NPROD')read(9,*)NPROD if(so(1:2).eq.'NW' )read(9,*)NW if(so(1:4).eq.'PRES')read(9,*)PRES if(so(1:4).eq.'TEMP')read(9,*)TEMP if(so(1:4).eq.'XMIN')read(9,*)xmin if(so(1:4).eq.'XMAX')read(9,*)xmax if(so(1:4).eq.'END')goto 99 goto 1 99 close(9) write(6,*)'Option read' write(6,*) inquire(file='ham.xyz.or',exist=lex) if(.not.lex)call report('ham.xyz.or not found!') inquire(file='ham.prm',exist=lex) if(.not.lex)call report('ham.prm not found!') inquire(file='KEY.TXT',exist=lex) if(.not.lex)call report('KEY.TXT not found!') c c List of options c --------------- c dynamic the dynamic command, with system path (e.g. /bin/dynamic) c minimize the minimize command, with system path c TINKER VERSION (ltinker5) c xmin lower coordinate limit c xmax upper coordinate limit c NEQUI number of equilibrium steps c NPROD number of production steps c NENS ensemble type (2-NVT, 4-NPT) c NFS MD integration step in fs c NPS MD time between dumps in ps --- only one dump per simulation c NINT histogram integration step -- not used here c NW number of weighted histograms c MTOL minimization gradient limit c PRES pressure (for NPT) c TEMP temperature c c KEY.TXT typically contains periodic boundary condition c definition (A-AXIS ...) and the restricted coordinate c open(9,file='goham') write(9,*)'cp ham.xyz.or ham.xyz' dx=(xmax-xmin)/dble(NW) x=xmin-dx write(9,*)'rm ham.xyz_* *.dyn ham.0* ham.1* ham.2* ham.3*' write(9,*)'rm ham.4* ham.5* ham.6* ham.7* ham.8* ham.9*' do 2 i=1,NW c write(so,'(i10)')i do 21 istart=1,10 21 if(so(istart:istart).ne.' ')goto 22 22 do 23 iend=10,1,-1 23 if(so(iend:iend).ne.' ')goto 24 24 x=x+dx c c prepare key file open(91,file='ham.key.'//so(istart:iend)) 3 open(92,file='KEY.TXT') read(92,980,end=33,err=33)s80 980 format(a80) c c periodic boundary conditions: if(s80(3:6).eq.'AXIS')lpc=.true. if(s80(1:17).eq.'RESTRAIN-DISTANCE')then read(s80(18:80),*)i1,i2,b write(91,9191)i1,i2,b,x,x 9191 format('RESTRAIN-DISTANCE',2i5,f16.8,2f12.4) else if(s80(1:16).eq.'RESTRAIN-TORSION')then read(s80(18:80),*)i1,i2,i3,i4,b write(91,9192)i1,i2,i3,i4,b,x,x 9192 format('RESTRAIN-TORSION',4i5,f16.8,2f11.2) else write(91,980)s80 endif endif goto 3 33 close(91) close(92) write(9,'(a)')'cp ham.key.'//so(istart:iend)//' ham.key' c c equilibration stage minimization and dynamics do 212 ist=1,80 212 if(minimize(ist:ist).ne.' ')goto 222 222 do 2332 ien=80,1,-1 2332 if(minimize(ien:ien).ne.' ')goto 2412 2412 write(9,9031)(minimize(j:j),j=ist,ien) 9031 format(1a1,$) write(9,9032)MTOL 9032 format(' ham ',f10.2) write(9,*)'mv ham.xyz_2 ham.xyz' do 211 ist=1,80 211 if(dynamic(ist:ist).ne.' ')goto 221 221 do 233 ien=80,1,-1 233 if(dynamic(ien:ien).ne.' ')goto 241 241 write(9,*)'rm CURRENT' write(9,9031)(dynamic(j:j),j=ist,ien) NPS=NFS*dble(NEQUI)/1000.0d0 if(lpc)then if(NENS.eq.4)then write(9,903)NEQUI,NFS,NPS,NENS,TEMP,PRES else write(9,903)NEQUI,NFS,NPS,NENS,TEMP 903 format(' ham ',i9,2f7.1,i2,2f7.1) endif else if(ltinker5)then write(9,9041)NEQUI,NFS,NPS,NENS,TEMP 9041 format(' ham ',i9,2f7.1,i2,f7.1) else write(9,904)NEQUI,NFS,NPS,TEMP 904 format(' ham ',i9,3f7.1) endif endif write(9,902) 902 format('mv ham.001 ham.xyz') c c production dynamics write(9,*)'echo ',i,' > CURRENT' write(9,9031)(dynamic(j:j),j=ist,ien) NPS=NFS*dble(NPROD)/1000.0d0 if(lpc)then if(NENS.eq.4)then write(9,903)NPROD,NFS,NPS,NENS,TEMP,PRES else write(9,903)NPROD,NFS,NPS,NENS,TEMP endif else if(ltinker5)then write(9,9041)NPROD,NFS,NPS,NENS,TEMP else write(9,904)NPROD,NFS,NPS,TEMP endif endif write(9,902) write(9,*)'cp ham.xyz '//so(istart:iend)//'.xyz' write(9,*)'rm *.dyn' 2 write(9,*) close(9) call report('goham written') end subroutine report(s) character*(*) s write(6,60)('*',i=1,len(s)) write(6,*) write(6,'(a)')s write(6,60)('*',i=1,len(s)) 60 format(1a1,$) write(6,*) stop end