tranal.f

分子动力学模拟程序

 	program TRANAL
*  Analys of trajectory file 
*  Calculations of orientational correlation functions 
*  three-body correlations and many other things
*
*  v.3.2    2005-11-22   Xmol trajectory: read box sizes from the 
*                        commentary line
*           2005-10-25   Correction of CH bond lengths
*  v 3.1    2005-04-22   Odred parameter /dipole coupling calculations
*  v.3.0    2005-02-10
*
	include "tranal.h"
        parameter (MAXFIL=1000)
	character*64 filnam(MAXFIL),fildif,filres,fildif2,
     +filtcf,fildfm,FFILTR,filxm,
     +PATHDB,FILPDB,FNAME,FANGL,FTOR,HOST*24,RNAME*96,TMPFIL*8,COM*128,
     +FXMOL,FIRDF,FRMS,CH*4,cdum*8
	dimension TIMIN(MAXFIL) 
	dimension OX(NTOT),OY(NTOT),OZ(NTOT)
	real*4 TX(NTOT),TY(NTOT),TZ(NTOT)
	logical logo,LORIENT,LCOR3,LDIFF,LRES,LPDB,LANG,
     +LRDF,LRMS,LTCF,LVEL,LFILTR,LORD
	namelist /INPUT/LCOR3,LORIENT,IS1,IS2,IS3,AMPL,fil3cor,filori
     +,RMAXX,NRX,RMAXY,NRY,ITH,ISOR,ROMAX,NOM,NOI,IO1,IO2,IO3,NFORM,
     +RB12,D12,IPRINT,NDUP,FSET,NOMX,NOMY,NOMZ,RXMAX,RYMAX,RZMAX,
     +ZB,DZB,IPROC,LDIFF,IDF,MTR,DTT,NTT,ISTEP,FILDIF,FILRES,ICNT,
     +LRES,IR1,IR2,RRN,RRX,RRD,NITT,TIMM,TTER,LCOL2,NDUR,NM,IDFM1,IDFM2,
     +LANG,FANGL,FTOR,NAA,LXMOL,FXMOL,RDFCUT,NA,NR1,NR2,LRDF,FIRDF,
     +FORDF,LPDB,PATHDB,NAMOL,NFBEG,NFEND,FNAME,FCRD,LORD,
     +FORD,FORDIN,SCCH,
     +LRMS,FRMS,IMAC,BOXL,BOYL,BOZL,NTYPES,NSPEC,IAT,ISELF,
     +LTCF,filtcf,N1,N2,ITCF,NSTEG,DTCF,LVEL,NAI,RMI,RMAX,
     +RB13,D13,RB23,D23,LFILTR,FFILTR,NUMTR,BREAKM
*  
*
*    Required parameters:
*
*    NFORM - trajectory format
*        MDYN     - binary MDynaMix
*        XMOL     - XMOL  (xyz)
*        PDBT     - PDB trajectory from GROMACS
*
*  LRDF - calculate many RDF-s
*    FIRDF - input file specifying rdf (the same type as in MDynamix)
*    FORDF - resulting rdfs
*    RDFCUT - RDF cutoff
*    NA    - RDF resolution
*    NAI   - RDF resolution for intramolek RDF
*    RMI,RMAX - min and max limits of intramolec RDF
*     
*  LORIENT - calculate orientational or spatial distribution functions
*    Recuired parameters:
*    NOI - number of equivalent local coord. systems 
*    IO1, IO2, IO3 (NOI) - define local coordinate system
*    IPROC = -1 only for DNA - SDF around two base pairs
*               IO1 = phosphates; IO2 - which counts; IO3 - which are
*               used to build average structure     
*             0 normal calculations of SDF defined by IO1,IO2,IO3
*             1,2,3  - slices SDF in 3 planes. may not work properly
*             4 only for DNA fragment : SDF around DNA as a whole
*               in this case DNA strands should be 1 and 2 type  
*    
*    FSET - file with atom numbers for calculation of average structure 
*
*  LDIFF - calculate diffusion
*    DTT - total time in s	
*    IAT - atom num. in the molecule for diffusion
*    ISTEP - interval
*    BREAKM (ps) - max break of trajectory 
*
*  LCOR3 - three-body correlations
*    IS1(NDUP),IS2(NDUP),IS3(NDUP) - site number
*    RMAXX, RMAXY - size of the output field
*    NRX,   NRY   - resolution
*    RB12  - fixed distance between 1 and 2
*    D12   - fluctuation of RB12
*
*  LANG - calculate distribution over torsion angles
*    Required parameters:
*    FANGL - name with atom numbers defining torsions
*    FTOR  - output file
*
*  LTCF - calculate tcf
*    N1(ITYP),N2(ITYP) - define vector on eachh molecule for reorient.tcf
*    ITCF(12) - which tcf calculate
*    NSTEG - num of steps for tcf
*    DTCF - time step of TCF
*
*  LORD - calculate order parameter / dipole couplings relative Z axis
*
*    FORDIN - file name containing input for order parameter calculations
*    FORD - file name containing output
*    SCCH - scale CH bonds by this factor
*
*  LRMS - root mean square deviation
*    FRMS - file for RMS
*    IMAC - for which molecule types calculate RMS (relative to .mol file) 
*    NMRMS - which molecule test for RMS (if 0 - all molecules of this type)
*    FILREF - filename with reference structure ( "XMOL" format)
*
*  other things:
*
*   LVEL - read velocities from the trajectory files
*   LFILTR - analysis for particular atoms
*   FFILTR - which atoms take for  the analysis
*   NUMTR - number of atoms in the trajectory
*
*  LXMOL - write trajectory in XMOL format
*    Required parameters:
*    FXMOL - name of files
*    ISTEP - how often dump comfigurations in rewrited trajectory
*
*  default values
	NFORM='XMOL'
	NDUP=1  
	NDUR=1
        ISTEP=1
	IAT=1
	NTT=200
	MTR=MTRM
	RRD=0.25
	RDFCUT=10.
	NA=100
	IS3(1)=1       
	IPRINT=5
	IPROC=0                     
	BREAKM=5.
	BOXL=1000.
	BOYL=1000.
	BOZL=1000.
	SCCH = 1.
	LORIENT=.false.
	LRDF=.false.
	LRES=.false.
	LPDB=.false.
	LANG=.false.
	LTCF=.false.
	LDIFF=.false.
	LFILTR=.false.
	LXMOL=.false.
	LVEL=.false.
        LORD=.false.
	do I=1,MAXCOLOR
	  AMPL(I)=0.2*I
	end do
	NRX=NMX
	NRY=NMY
	ITH=1
	ISOR(1)=1
	ROMAX=5.
	NOM=50
	NOMX=50
	NOMY=50
	NOMZ=50
	ZB=0.
	DZB=0.1 
	NAA=36
	ICTCF=0
	fil3cor='correl3.ps'
	filori='orient.ps'
	fildif='tcor.dat' 
        filres='restime.lst' 
	fcrd='mean.crd'
	filpdb='fort.7' 
	ftor='tordist.lst'
	fxmol='fxmol.lst'
	frms='frms.dat'
* Data input
	J=0             
	NTP=0 
	IADT(1)=0   
	read(*,INPUT)                              
	NF=NFEND-NFBEG+1       
	if(NFORM.eq.'PDBT')LVEL=.false.
* Parameters for analysis of dihedrals  
	if(LANG)then
	  open(unit=17,file=FANGL,status='old',err=199)
	  IE=0 
	  STR=TAKESTR(17,IE)
	  read(STR,*)NAG               ! num. of types 
	  if(NAG.gt.NTOR)stop ' increase NTT'
	  if(NA.gt.NAAM)stop ' increase NAAM'
	  STR=TAKESTR(17,IE)
	  read(STR,*)(MANGL(I),I=1,NAG)    ! torsions of each type
	  IADT(1)=0
	  do I=1,NAG
	    IADT(I+1)=IADT(I)+MANGL(I)
	  end do
	  NADT=IADT(NAG+1)                   ! total  
	  if(NADT.gt.NT)stop 'increase NT'
	  do I=1,NADT                   
	    STR=TAKESTR(17,IE)
  	    READ(STR,*)IT(I),JT(I),KT(I),LT(I)    ! which
	  end do
	  close(17)
	  write(*,*)' atoms for torsion distributions are read' 
	end if
	NCF=0
	if(LORIENT.and.(IPROC.le.0.or.IPROC.ge.4))then
*   coordinates of which sites will be averaged in local coord.system
	   write(*,*)NOMZ,NOMY,NOMX
	   write(*,'(6e13.5)')RXMAX,RYMAX,RZMAX
	   do II=1,NOI
	     if(IO3(II).gt.0)then
		NCF=NCF+1
		ICF(II)=NCF
	     end if
	   end do
	   write(*,*)NCF,(ICF(I),I=1,NCF)
	   if(NCF.eq.0.and.IPROC.le.0)go to 8
	   open(unit=14,file=FSET,status='old',err=901)
	   read(14,*)NMSTR
	   do I=1,NMSTR
	      read(14,*,err=902,end=902)ISREP(I),(ISM(I,II),II=1,NCF)
	   end do
	   close(14)
	   write(*,*)' read ',NMSTR,' reference atom coordinates'
	   go to 8
 901	   write(*,*)' file ',FSET,' not found'
	   stop
 902	   write(*,*)' input error in file ',FSET,' line ',I
	   stop
	end if
 8	if(LXMOL)then
 	  LLNX=0 
 	  do I=1,64
 	    ISTR=ichar(FXMOL(I:I))
	    if(ISTR.le.15)then
	      LLNX=I-1
	      go to 18
	    end if
            if(FXMOL(I:I).ne.' ')LLNX=I
	  end do 
	end if
*    check array dimensions
 18	if(NAA.gt.NTA)stop ' increase NTA'
	if(NF.le.0)stop '!!! invalid NFBEG,NFEND'  
	if(NFEND.ge.1000)stop '!!! NFEND > 999 not supported' 
	if(NF.gt.MAXFIL)stop '!!! increase MAXFIL !!!'
	if(NOI.gt.NDUPM)stop '!!! increase NDUPM'
	if(NOMX.gt.NOMXM)stop 'increase NOMXM !!!'
	if(NOMY.gt.NOMYM)stop 'increase NOMYM !!!'
	if(NOMZ.gt.NOMZM)stop 'increase NOMZM !!!'
	if(NA.gt.NAAM)stop 'increase NAAM !!!'
	IL=0 
	do I=1,64
	  ISTR=ichar(FNAME(I:I))
	  if(ISTR.le.15)then
	    LF=I-1
	    go to 11
	  end if
        if(FNAME(I:I).ne.' ')IL=I
	end do 
	LF=IL
*  Check and reconstruct file list
*  Check that the file list is on a remote computer
 11	JJ=0
	do I=1,LF-1
	  if(FNAME(I:I).eq.'/')JJ=I
	end do	
 26	LFJ=LF-JJ
	write(*,*)FNAME(1:LF)
	J=0
 	do I=1,NF
	  NFL=NFBEG+I-1
	  filnam(I)(1:LF)=FNAME(1:LF)
	  filnam(I)(LF+1:LF+4)='.000'
	  do M=LF+5,64
	    filnam(I)(M:M)=' '
	  end do
	  if(NFL.lt.10)then
	    write(filnam(I)(LF+4:LF+4),'(i1)')NFL
	  else if(NFL.lt.100)then
	    write(filnam(I)(LF+3:LF+4),'(i2)')NFL
	  else
	    write(filnam(I)(LF+2:LF+4),'(i3)')NFL
	  end if 
	  write(*,*)' looking for file ',filnam(I)(1:40)
	  if(NFORM.eq.'PDBT')then
            open (unit=12,file=filnam(I),status='old',err=5)
	if(IPRINT.ge.5)write(*,*)' file ',filnam(I),' is open'
 9	    read(12,'(a)',err=10)STR
*     PDB trajectory format
	    if(STR(1:6).ne.'HEADER'.and.STR(1:5).ne.'TITLE')go to 9
	    do K=7,80
	      if(STR(K:K+1).eq.'t=')then
	        read(str(K+2:K+20),*)FULTIM
  	        if(IPRINT.ge.5)write(*,*)' FULTIM ',FULTIM,' ps' 
		go to 13
	      end if
	    end do
	    write(*,*)' time is not defined '
 13	    continue
	    go to 20 
*  XMOL trajectory format
	  else if(NFORM.eq.'XMOL')then
            open (unit=12,file=filnam(I),status='old',err=5)
*	if(IPRINT.ge.5)write(*,*)' file ',filnam(I),' is open'
	    read(12,*,err=10)NSTOT
	    read(12,*,err=10)cdum,FULTIM
	    FULTIM=0.001*FULTIM            ! in ps
	    if(IPRINT.ge.5)write(*,*)' FULTIM ',FULTIM,' ps'  
 	    go to 20 
*  Binary MDynaMix format
	  else if(NFORM.eq.'MDYN')then
            open (unit=12,file=filnam(I),status='old',
     +    form='unformatted',err=5)
	    if(IPRINT.ge.5)write(*,*)' file ',filnam(I),' is open'
            read(12,err=10) IA, dt, boxl, boyl, bozl, unitl, ntypes
	    if(IA.eq.1)LVEL=.true.
	    if(IPRINT.ge.5)write(*,*)' num.of types ',NTYPES  
	    if(I.gt.1)then
	      if(NTP.ne.NTYPES)then
	        write(*,*)
     +      ' uncompatible NTYPES in ',filnam(I),' and ',filnam(I-1) 
	        stop
	      end if
	    else
	      rewind (12)
              read(12,err=10) IA, dt, boxl, boyl, bozl, unitl, ntypes,
     x                (nspec(ii),nsits(ii),logo,ii=1,ntypes)
	    end if 
	    NTP=NTYPES  
	    do ITYP=1,NTYPES
	      read(12,err=10)
	      if(IPRINT.ge.8)write(*,*)ITYP,NSPEC(ITYP),NSITS(ITYP)
	    end do
	    read(12)IA,FULTIM
	    FULTIM=FULTIM*1.d12         ! in ps
	    if(IPRINT.ge.5)write(*,*)' FULTIM ',FULTIM,' ps'  
	    close(12)
	    go to 20 
	  else
	    write(*,*)' Unsupported trajectory format ',NFORM
	    stop
	  end if
  5	  write(*,*)' File not found: ',filnam(I)
	  go to 20
 10	  write(*,*)' File ',filnam(I),': input error'
 20	  close(12)
	  J=J+1
	  TIMIN(J)=FULTIM 
	  filnam(J)=filnam(I)
	end do
      write(*,*)' input of the file list is completed'	  
 30	NF=J
	FULTIM0=TIMIN(1)
 	TTER=TTER*1.d-12    
	if(LORIENT)then
	  if(NOMAX.gt.NMX)stop 'too large NOMAX'
	  if(NOMAX.gt.NMY)stop 'too large NOMAX'
	  if(NOM.gt.NOMAX)stop 'increase NOMAX' 
	  if(IPROC.eq.0.and.2*NOM.gt.NOMAX)stop 'increase NOMAX'
	end if 
	if(LDIFF)then
	  if(MTR.gt.MTRM)stop ' increase MTRM'
	  if(NTT.gt.NTCD)stop ' increase NTCD'
	  do ITYP=1,NTT
	    RCUM(ITYP)=0.d0
	    NCUM(ITYP)=0
	  end do   
	  do I=1,MAXMOL
	    ISHX(I)=0
	    ISHY(I)=0
	    ISHZ(I)=0
	  end do
	end if 
	if(LRES)then  
	  if(NITT.gt.NITTM)stop 'increase NITTM'
	  if(RRN.gt.RRX)stop 'RRN>RRX'
	  RRNN=RRN-RRD
	  RRNX=RRN+RRD
	  RRXN=RRX-RRD
	  RRXX=RRX+RRD
	end if
	if(LRMS)then
	   open(unit=24,file=frms,status='unknown')
	end if
	if(NRX.gt.NMX)stop ' increase NMX' 
	if(NRY.gt.NMY)stop ' increase NMY'
*  for 2D distributions 
	RMX12=(RB12+0.5*D12)**2
	RMN12=(RB12-0.5*D12)**2 
	DXX=2.d0*RMAXX/NRX
	DYY=2.d0*RMAXY/NRY
	DRI=RMAX-RMI
*  for 3D distributions
	DDX=2.d0*RXMAX/NOMX
	DDY=2.d0*RYMAX/NOMY
	DDZ=2.d0*RZMAX/NOMZ
	ROMAX2=(RXMAX**2+RYMAX**2+RZMAX**2)
	DDO=ROMAX/NOM     
	RDF2X=(RDIF2+DDIF2)**2
	RDF2N=(RDIF2-DDIF2)**2	
*  Input other parameters
*  Adresses and referenses
	IAN=0 
	ISADR(1)=0
	ISADDR(1)=0    
	IADDR(1)=0
	NOP=0 
	NSITES=0
	NSTOT=0 
	FNSTR=0.
	IATOM=0
	IMOL=0
	IST=0   
	NCOR=0    
	JBREAK=1                      
	NOM2=2*NOM
	RMSS = 0.
	TRTEMP=0.
	NX=0
	do ITYP=1,NTYPES
	  LIST(ITYP)=1
	  if(NFORM.eq.'MDYN')NSTS=NSITS(ITYP)
          CALL READMOL (SUMM,NX,ITYP,PATHDB,
     +  NAMOL(ITYP))
	  if(NFORM.eq.'MDYN'.and.NSTS.ne.NSITS(ITYP))then
	    write(*,*)' Molecule of type ',ITYP,' has ',NSTS,
     &' sites in the trajectory but ',NSITS(ITYP),' in the input file'
	    stop
	  end if
          TOTMAS       = TOTMAS+SUMM*NSPEC(ITYP)
          SUMMAS(ITYP) = SUMM 
	  NAME(ITYP)=NAMOL(ITYP)(1:6)
        end do                
        FACTM=AVCNO*1.d3
        do ityp=1, ntypes
	  ISADR(ITYP+1)  = ISADR(ITYP)+ NSITS(ITYP)          ! site addr
	  ISADDR(ITYP+1) = ISADDR(ITYP)+NSPEC(ITYP)*NSITS(ITYP)  ! atom adr
          IADDR (ITYP+1) = IADDR (ITYP)+NSPEC(ITYP)         ! mol. addr
          NOP     = NOP    + NSPEC(ITYP)
          NSITES  = NSITES + NSITS(ITYP)
          if(NSITS(ITYP).eq.2)FNSTR=FNSTR+2.d0*NSPEC(ITYP)
          if(NSITS(ITYP).ge.3)FNSTR=FNSTR+3.d0*NSPEC(ITYP)
          NSTOT   = NSTOT  + NSPEC(ITYP)*NSITS(ITYP)
          FNSS3   = 3.D0*DFLOAT(NSITES)
          write(*,'(a6,i4,3x,a5,i5,4x,a6,i4)')
     +' type ',ITYP,' Nmol',NSPEC(ITYP),
     +'Nsites',NSITS(ITYP)
	  do J=ISADR(ITYP)+1,ISADR(ITYP+1)
	    MASSD(J)=MASS(J)/TOTMAS/FACTM
	    ITS(J)=ITYP
	  end do                     !    Nsite -> type
	  do IS=1,NSPEC(ITYP)
	    IMOL=IMOL+1
            do I=1,NSITS(ITYP)
	      ISITE=ISADR(ITYP)+I
	      IATOM=IATOM+1
	      NNUM(IATOM)=IMOL        !   Natom -> Nmol
	      NSITE(IATOM)=ISITE      !   Natom -> Nsite
	      ITYPE(IATOM)=ITYP       !   Natom -> type
	      if(IPRINT.ge.9)write(*,'(4I6,2x,2a4,a6)')
     +        IATOM,IMOL,ISITE,ITYP,NM(ISITE)
	    end do
	  end do
        end do                   
* Renumeration array
	if(LFILTR)then