tranal.f

分子动力学模拟程序

	   do I=1,NUMTR
	     NUMR(I)=0
	   end do
	   write(*,*)' renumerating atoms from file ',FFILTR
	   open(unit=14,file=FFILTR,status='old')
	   do I=1,NSTOT
	      read(14,*,err=921,end=921)II,NUMO(II)
	      NUMR(NUMO(II))=II
	      if(IPRINT.ge.7)write(*,*)I,II,NUMO(II)
	   end do
	   close(14)
	else
	   NUMTR=NSTOT
	   do I=1,NSTOT
	      NUMO(I)=I
	      NUMR(I)=I
	   end do	   
	end if
*  number of torsions to compute
	if(LANG)then
	  do I=1,NAG
	    MANGL(I)=MANGL(I)*NSPEC(ITS(IT(IADT(I)+1)))
	  end do
	end if
*   Input RDFs
	BOXS=0.           ! to accumulate average box sizes
	BOYS=0.
	BOZS=0.
	if(LRDF)then
	  open(unit=17,file=FIRDF,status='old',err=197)
	  IE=0
	  N=0
	  STR=TAKESTR(17,IE)
	  read(STR,*)NRDF
	  write(*,*)' Number of RDF-s: ',NRDF
	  if(NRDF.gt.MAXRDF)stop 'Increase MAXRDF!!!'
          DO K       = 1,NRDF
 	    IE=-1
	    STR=TAKESTR(17,IE)
	    if(str(1:1).eq.'&')then
              READ(str(2:24),*)IRP
	      if(IRP.gt.MAXGR) stop 'Increase MAXGR!!!'
	      IREP(K)=IRP
	      do IR=1,IRP
	        STR=TAKESTR(17,IE)
	        read(STR,*,err=99)I,J  
		IRD(K,IR)=I
		JRD(K,IR)=J
	      end do
	    else 
	      IREP(K)=1
	      read(STR,*,err=99)I,J  
	      IRD(K,1)=I
	      JRD(K,1)=J
	    end if
	    write(*,*)' RDF  No ',K,' equiv ',IREP(K)
	    do I=1,IREP(K)
	      write(*,'(2I6,2(3x,a4))')
     +        IRD(K,I),JRD(K,I),NM(IRD(K,I)),NM(JRD(K,I))
	    end do
          END DO
	  go to 101
 99	  write(*,*)' input error in file ',FIRDF
	  stop
 101	  continue
	end if
*  Compose list of atoms for calculation of residence time
	if(LRES)then
	   NRT1=0
	   NRT2=0
	   IHP=0
	   do I=1,NR1
	      IS=IR1(I)
	      ITYP1=ITS(IS)
	      NS1=NSPEC(ITYP1)
	      do J=1,NS1
		 NRT1=NRT1+1
		 if(NRT1.gt.NRT1M)stop ' increase IRT1M'
		 ISP=ISADDR(ITYP1)+(J-1)*NSITS(ITYP1)+IS-ISADR(ITYP1)
		 IRS1(NRT1)=ISP
	      end do
	   end do
	   do I=1,NR2
	      IS=IR2(I)
	      ITYP2=ITS(IS)
	      NS2=NSPEC(ITYP2)
	      do J=1,NS2
		 NRT2=NRT2+1
		 if(NRT2.gt.NRT2M)stop ' increase IRT2M'
		 ISP=ISADDR(ITYP2)+(J-1)*NSITS(ITYP2)+IS-ISADR(ITYP2)
		 IRS2(NRT2)=ISP
	      end do
	   end do
	   if(IPRINT.ge.8)write(*,*)NRT1,(IRS1(I),I=1,NRT1)
	   if(IPRINT.ge.8)write(*,*)NRT2,(IRS2(I),I=1,NRT2)
	end if	
* check list of torsions
	if(IPRINT.ge.6.and.LANG)then 
	  write(*,*)' distribution over torsion angles:'
	  do I=1,NAG
          write(*,*)I,' type'  
	    do J=IADT(I)+1,IADT(I+1)  
	      write(*,*)NM(NSITE(IT(J))),NM(NSITE(JT(J))),
     + NM(NSITE(KT(J))),NM(NSITE(LT(J)))
	    end do
	  end do               
	end if
	if(NFORM.eq.'XMOL')then
	  HBOXL=0.5*BOXL
	  HBOYL=0.5*BOYL
	  HBOZL=0.5*BOZL
	  VOL=BOXL*BOYL*BOZL 
	end if
*  Zero tables
	do I=1,NA
	  RDF(I)=0.
	end do
	do IX=0,NRX
	  do IY=0,NRY
	    P3COR(IX,IY)=0.d0
	    P3CO2(IX,IY)=0.d0
	  end do
	end do
	do IX=0,2*NOM
	  do IY=0,NOM
	    POR(IX,IY)=0.
	  end do
	end do
	do I=1,NAG
	  IGS(I)=0
	  ITRAN(I)=0
	  do J=1,NAA
	    NPOP(I,J)=0 
	  end do
	end do 
	do IX=0,NOMXM
	   do IY=0,NOMYM
	      do IZ=0,NOMZM
		 ID3(IX,IY,IZ)=0
	      end do
	   end do
	end do
	do I=1,NS
	   XAV(I)=0.
	   YAV(I)=0.
	   ZAV(I)=0.
	   XAV2(I)=0.
	   YAV2(I)=0.
	   ZAV2(I)=0.
	   IAA(I)=0
	end do
	ITIR=0
	ITIR1=0
	ITIR2=0
	TIMRS=0
	TIMRS1=0.
	TIMRS2=0.    
	IPDB=6
	COR2=0.d0
	COR22=0.d0                 
	FULTIM1=TIMIN(1)
	write(*,*)' FULTIM1=',FULTIM1
	if(LDIFF)write(*,*)
     +' Calculation of the diffusion coefficients for type ',IDF,
     + NSPEC(IDF),NSITS(IDF)
*  Initialize TCF
	if(LTCF)call TCFINP
*  Initialize order parameter calculations
	if(LORD)call ORDERIN
*  Put other initialization things here:
*
*  Begin reading of next trajectiry file 
	ITRI=0
	do IFL=1,NF        
	  IIN=0     ! Signals that we read headers
*  Reading HEADERs
*  Case of PDB trajactory
 35	  if(NFORM.eq.'PDBT')then
	    if(IIN.eq.0)then
              open (unit=12,file=filnam(IFL),status='old',
     +        form='formatted') 
	      if(IPRINT.ge.7)write(*,*)' file ',IFL,' open'
	    end if
 31	    read(12,'(a)',end=45) STR
	    if(STR(1:6).ne.'HEADER'.and.STR(1:5).ne.'TITLE')go to 31 
	    do K=7,80
	      if(STR(K:K+1).eq.'t=')then
	        read(str(K+2:K+20),*)FULTIM
		go to 33
	      end if
	    end do
	    write(*,*)' time is not defined '
 33	    IFLR=0
* Line 2 
 32	    read( 12,'(a)',end=45 ) STR
	    if(STR(1:4).eq.'ATOM')then
	      if(IPRINT.ge.6)write(*,*)' box sizes undetermined'
	      IFLR=1
	      go to 34
	    end if
	    if(STR(1:4).ne.'CRYS')go to 32
	    read(str,*)cdum,BOXL,BOYL,BOZL
	    if(IPRINT.ge.6)write(*,'(a,3f12.3)')
     &      ' Box: ',BOXL,BOYL,BOZL
	  else if(NFORM.eq.'XMOL')then
*  XMOL trajectory
	    if(IIN.eq.0)then
              open (unit=12,file=filnam(IFL),status='old',
     +    form='formatted') 
	      if(IPRINT.ge.7)write(*,*)' file ',IFL,' open'
	    end if
 	    read(12,*,err=40,end=45) NUMTR
            read( 12,'(a)',err=40,end=40 )STR
            read(STR,*,end=40,err=40)cdum,FULTIM
	    do J=1,80
              if(STR(J:J+2).eq.'BOX')then
                read(STR(J+4:80),*,end=40,err=40)BOXL,BOYL,BOZL
	        if(IPRINT.ge.6)write(*,'(a,3f12.3)')
     &          ' Box: ',BOXL,BOYL,BOZL
	        go to 37
	      end if
	    end do
	    if(IPRINT.ge.7)write(*,*)' box sizes undetermined'
 37	    FULTIM=0.001*FULTIM	! in ps
	  else
*  Binary MDynamix trajectory
	    if(IIN.eq.0)then
              open (unit=12,file=filnam(IFL),status='old',
     +        form='unformatted') 
              read(12) IVEL, dt, boxl, boyl, bozl, unitl, ntypes,
     x                (nspec(i),nsits(i),logo,i=1,ntypes)
	      if(LXMOL)then
	        filxm(1:LLNX)=FXMOL(1:LLNX)
	        filxm(LLNX+1:LLNX+4)=filnam(IFL)(LF+1:LF+4)
	        open(unit=19,file=filxm,status='unknown')
	      end if
	      if(IPRINT.ge.7)write(*,*)
     &          'read first line.',NTYPES,' types'
              do ityp=1, ntypes
                NSB = ISADR (ITYP)+1
                NSE = ISADR (ITYP+1)
* Lines from 2 to NTYPES+1
                read( 12 ) (mass(j), j=NSB,NSE),LIST(ITYP)
		if(IPRINT.ge.8)write(*,*)ityp,NSB,NSE,LIST(ITYP)
              end do                    
	      if(IPRINT.ge.7)write(*,*)' read masses   '
	    end if
* First line of a configureation
            read(12,err=40,end=45)IA,fultim,cumtim,
     &      TEMP,PRES,EP,BOXL,BOYL,BOZL,(LIST(LL),LL=1,NTYPES)
	    FULTIM=FULTIM*1.d12                 ! in ps
	    if(IPRINT.ge.6)write(*,'(2(a,3f10.2))')
     &      'Box: ',BOXL,BOYL,BOZL,' term:',TEMP,PRES,EP
	  end if    ! if(NFORM...
 34	  if(IIN.eq.0)then
 	    write(*,'(2a/a,f14.3,a8,6I3)')
     +  ' begin analys of file  ',filnam(IFL), 
     +  ' init. T = ',fultim
	    IIN=1
	  end if  
	  IBREAK=0
	  DIFT=FULTIM-FULTIM1
  	  HBOXL=0.5*BOXL
	  HBOYL=0.5*BOYL
	  HBOZL=0.5*BOZL
	  VOL=BOXL*BOYL*BOZL 
	  if(DIFT.gt.BREAKM)then
	    write(*,*)
     +    'break in trajectory file. from ',FULTIM1,' to ',
     +    FULTIM,'ps'
	    IBREAK=1
	  end if    
	  if(IPRINT.ge.5.and.(LCOR3.or.LORIENT))
     +    write(*,'(a,f10.3,3x,a,f12.0,I12)')
     +    ' time=',fultim,' N conf.',COR2,NCOR
	  if(IFL.lt.NF.and.fultim.ge.TIMIN(IFL+1))go to 45
	  FULTIM1=FULTIM 
	  if(IPRINT.ge.7)write(*,*)' proceeding time ',FULTIM
	  if(NFORM.eq.'MDYN')then
c - for each type
            do ityp=1,ntypes
              NSB          = ISADDR(ITYP)+1
              NSE         = ISADDR(ITYP+1) 
* Other NTYPES lines of each point
              if(LIST(ITYP).ne.0)then
	        if(IVEL.le.0)then
	          read(12,err=40,end=45) (TX(j),j=NSB,NSE),
     *                (TY(j),j=NSB,NSE),
     *                (TZ(j),j=NSB,NSE)           
 	          do J=NSB,NSE
	            SX(J)=TX(J)
	            SY(J)=TY(J)
	            SZ(J)=TZ(J)
	          end do
	        else    !  IVEL=1
	          read(12,err=40,end=45) (TX(j),j=NSB,NSE),
     *                (TY(j),j=NSB,NSE),
     *                (TZ(j),j=NSB,NSE)
 	          do J=NSB,NSE
	            SX(J)=TX(J)
	            SY(J)=TY(J)
	            SZ(J)=TZ(J)
	           end do
	           read(12,err=40,end=45) (TX(j),j=NSB,NSE),
     *                (TY(j),j=NSB,NSE),
     *                (TZ(j),j=NSB,NSE)
 	           do J=NSB,NSE
	             VX(J)=TX(J)
	             VY(J)=TY(J)
	             VZ(J)=TZ(J)
	           end do
	         end if    !  IVEL
	       end if     ! LIST
             end do    !  ITYP
	  else     ! ASCII trajectories
            do i=1,NUMTR
	      OX(I)=SX(I)
	      OY(I)=SY(I)
	      OZ(I)=SZ(I)
	      IN=NUMR(I)
	      if(LVEL)then
	        read(12,*,err=40,end=40)CH,XX,YY,ZZ,VVX,VVY,VVZ
		if(IN.ne.0)then
		   SX(IN)=XX
		   SY(IN)=YY
		   SZ(IN)=ZZ
		   VX(IN)=VVX
		   VY(IN)=VVY
		   VZ(IN)=VVZ
		   if(IPRINT.ge.8)write(*,'(2i5,3f10.4)')IN,I,XX,YY,ZZ
		end if
              else
		if(NFORM.eq.'PDBT')then
 53		  read(12,'(a)',end=40)STR
		  if(STR(1:6).eq.'ENDMDL'.or.STR(1:6).eq.'HEADER')then
		    write(*,*)' Number of atoms in the trajectory ',II,
     &		    ' < than in the input file ',NSTOT
		    stop
		  end if
		  if(STR(1:4).ne.'ATOM')go to 53
		  read(str,*,err=40)cdum,II,CH,cdum
		  read(str(31:80),*)XX,YY,ZZ
		else
	          read(12,*,err=40,end=40)CH,XX,YY,ZZ
		end if		
		if(IN.ne.0)then
		   SX(IN)=XX
		   SY(IN)=YY
		   SZ(IN)=ZZ
		   if(IPRINT.ge.8)write(*,'(2i5,3f10.4)')IN,I,XX,YY,ZZ
		end if
	      end if
	  end do
	end if
        if(SCCH.lt.0.9999.or.SCCH.gt.1.00001)call CHSCALE
	JBREAK=0
	if(ITRI.ne.0.and.DIFT.le.0.00001)then
	   write(*,*)' double point in the trajectory file at ',FULTIM
	   go to 35
	end if
	ITRI=1
	DT=DIFT
	if(DIFT.le.0.0001)DIFT=0.0001
	EKIN=0.
	do I=1,NSTOT
	  DX=SX(I)-OX(I)
	  DY=SY(I)-OY(I)
	  DZ=SZ(I)-OZ(I)
 	  if(DX.gt.HBOXL)DX=DX-BOXL
 	  if(DX.lt.-HBOXL)DX=DX+BOXL
 	  if(DY.gt.HBOYL)DY=DY-BOZL
 	  if(DY.lt.-HBOYL)DY=DY+BOZL
 	  if(DZ.gt.HBOZL)DZ=DZ-BOZL
 	  if(DZ.lt.-HBOZL)DZ=DZ+BOZL
	  IS=NSITE(I)
	  if(.not.LVEL)then
	    VX(I)=DX/DIFT
	    VY(I)=DY/DIFT
	    VZ(I)=DZ/DIFT
	  end if
	  EKIN=EKIN+MASS(IS)*(VX(I)**2+VY(I)**2+VZ(I)**2)   ! 2 Ekin
	  OX(I)=SX(I)
	  OY(I)=SY(I)
	  OZ(I)=SZ(I)
	end do
	TEMP=EKIN/(0.3*AVSNO*NSTOT*BOLTZ)
*
*  At this point, arrays SX,SY,SZ contains current coordinates of atoms
*                        VX,VY,VZ - their velocities
*
*  Calculations of different properties begin
*
	  if(LCOR3)then
	    do IGR=1,NDUP   
	      if(IGR.le.NDUR)then
	        call COR3(IGR,1)
	      else
		call COR3(IGR,2)
	      end if
	    end do
	  end if
	  if(LORIENT)then
	    do IGR=1,NDUP
	      call ORIENT(IGR)
	    end do
	  end if   
	  if(LANG)then
	    call  TAKETOR
	  end if     
	  if(mod(IAN/NDUP,ISTEP).eq.0.and.LDIFF)then
	    call DIFFUS
	  end if             
	  if(LPDB.and.mod(IAN/NDUP,ISTEP).eq.0)then
	    IPDB=IPDB+1
	    call PDBDUMP(IPDB)
	  end if
	  if(LRES)then
	    call RESTIM(0)
	  end if
	  if(LRDF)then
	     call RDFCALC
	  end if
	  if(LTCF.and.(FULTIM-FULTIM0)+1.d-6.ge.DTCF*ICTCF)then
 237	     ICTCF=ICTCF+1
	     call GETCOM
	     call DIPOLE
	     call GETROT
	     call GETTCF
	     TRTEMP=TRTEMP+TEMP
	     if(IPRINT.ge.6)write(*,'(a,I6,f11.5,3x,f10.2)')
     +  ' tcf: ',ICTCF,FULTIM,TEMP
	     if(FULTIM-FULTIM0.gt.DTCF*ICTCF)then
		write(*,*)' no config. at t= ',DTCF*ICTCF+FULTIM0
		go to 237
	     end if
	     if(IPRINT.ge.8)then
	       write(*,*)' local atom coordinates in principal s'
	       do I=1,NSTOT
		 IMOL=NNUM(I)
		 QX(IMOL)=WX(I)
		 QY(IMOL)=WY(I)
		 QZ(IMOL)=WZ(I)
		 call ROTMAT(QX,QY,QZ,IMOL,1,.TRUE.)
             write(*,'(2I5,3f10.4)')I,IMOL,QX(IMOL),QY(IMOL),QZ(IMOL)
	       end do
	     end if
	  end if
	  if(mod(IAN/NDUP,ISTEP).eq.0.and.LRMS)then
	     call RMSCALC
	  end if
	  if(LORD)then
	     call ORDERP
	  end if
*  Put your extra subroutines here
*
*         call SOMETHING
*
*  Dump trajectory of selected atoms in XMOL format, with ISTEP
	  if(mod(IAN/NDUP,ISTEP).eq.0.and.LXMOL)then
	     NSAT=0.