tranal.f

分子动力学模拟程序

	     do ITYP=1,NTYPES
	       if(LIST(ITYP).ne.0)NSAT=NSAT+NSITS(ITYP)*NSPEC(ITYP)
	     end do
	     write(19,*)NSAT
	     write(19,*)' after ',FULTIM*1.d3,' fs'
	     do I=1,NSTOT
		ITYP=ITYPE(I)
		IS=NSITE(I)
		if(LIST(ITYP).ne.0)write(19,'(a2,3(1x,f9.4))')
     +          NM(IS)(1:2),SX(I),SY(I),SZ(I)
	     end do
	  end if
	  IAN=IAN+NDUP
	  go to 35
 40	  write(*,*)' input error in file ',filnam(IFL),' T=',FULTIM
 45	  write(*,'(2a/a,f10.3,a,f12.0,i10,2i8)')
     +' finish analys of file ',filnam(IFL), 
     +  ' final T = ',fultim,'      N conf.  ',COR2,NCOR,ITIR,NDP2
	  close(12)
	end do
	write(*,*)' finish analys of the trajectory'
      write(*,*)' total number of configurations ',IAN  
*   Base pair RDF
	if(LCOR3)then  
	  NSS=NSPEC(ITS(IS1(1)))
	  NPAIRS=NSS*NSPEC(ITS(IS2(1)))
	  if(IS1(1).eq.IS2(1))NPAIRS=NSS*(NSS-1)
	  SH12=4*PI*D12*(RB12**2+D12**2/12.d0)
	  RDF2=COR2*VOL/(dfloat(IAN*NPAIRS)*SH12)
	  write(*,'(a,f8.3,a,f9.5)')
     +' 1-2 pair correlation at r=',RB12,' is ',RDF2
	end if
*   Final output
	if(LCOR3)call MCOR3
	if(LORIENT)call MORI
	if(LDIFF)call MDIFF(fildif)
	if(LRDF)call RDFOUT
	if(LORD)call ORDEROUT
	if(LTCF)then
	  TRTEMP=TRTEMP/ICTCF
	  call TCFOUT(filtcf)
	end if
	if(LRES)then
	  call RESTIM(1)
          open(unit=15,file=filres,status='unknown')
	  write(15,*)'#  distribution over time residence '
	  write(15,*)'#  sites ',IR1(1),' and ',IR2(1)
	  write(15,*)'#  between ',RRN,' and ',RRX,' +/- ',RRD 
	  write(15,*)'#  total ',ITIR2,' and ',ITIR,' events'
	  write(15,*)'#  cases of longer T* time ',ITIR1
	  write(15,*)'#  permission time T* = ',TTER*1.d12,' ps '
	  write(15,*)'#  average number of hydrated pairs ',
     +           IHP*1.*NDUP/IAN
	  write(15,*)
	  write(15,*)'#  av. residence time with T* ',
     +          TIMRS2/ITIR2,' ps' 
	  write(15,*)'#  simple average residence time',TIMRS/ITIR,' ps' 
	  write(15,*)'#  exponent decay            ',
     +  TIMRS1/ITIR1-TTER*1.d12,' ps'
	  write(15,*)
	  write(15,*)'#  time(ps)      hits1      P1     DECAY    ',
     +  ' hits2      P2'
	  DTIRS=ITIR2
	  do I=1,NITT  
	    PPP1=DIST2(I)/ITIR2
	    DTIRS=DTIRS-DIST2(I)
	    PFUN=DTIRS/ITIR2          
	    TTT=I*TIMM/NITT                 
	    PPP=DISTT(I)/ITIR          
	    write(15,'(f12.3,2x,f6.0,2f13.8,2x,f6.0,f13.8)')
     +    TTT,DIST2(I),PPP1,PFUN,DISTT(I),PPP
	  end do
	end if  
	if(LANG)then
	  open(unit=17,file=ftor,status='unknown')
	  write(17,'(a)')'# distribution over torsion angles'
	  write(17,'(a1,10x,50(a4,5x))')'#',
     +(NM(NSITE(IT(IADT(I)+1))),I=1,NAG)
	  write(17,'(a1,10x,50(a4,5x))')'#',
     +(NM(NSITE(JT(IADT(I)+1))),I=1,NAG)
	  write(17,'(a1,10x,50(a4,5x))')'#',
     +(NM(NSITE(KT(IADT(I)+1))),I=1,NAG)
	  write(17,'(a1,10x,50(a4,5x))')'#',
     +(NM(NSITE(LT(IADT(I)+1))),I=1,NAG)
	  do M=1,NAA
	    ANG=-180+(M-0.5)*360/NAA 
	    FACT=NAA*NDUP*1.d0/IAN
	write(17,'(f7.1,50f9.4)')ANG,
     + (NPOP(I,M)*FACT/MANGL(I),I=1,NAG)
	  end do
	  write(17,'(a,50f9.6)')'# Gaushe ',
     +(IGS(I)*1./(IGS(I)+ITRAN(I)),I=1,NAG)
	  write(17,'(a,50f9.6)')'# Trans  ',
     +(ITRAN(I)*1./(IGS(I)+ITRAN(I)),I=1,NAG)
	  close(17)
	end if
	if(LRMS)then
	   write(*,*)' R M S = ',sqrt((RMSS*ISTEP)/(IAN*NDUP))
	   write(24,'(a,f14.5)')'#    R M S = ',
     +     sqrt((RMSS*ISTEP)/(IAN*NDUP))
	end if
	stop
 197	write(*,*)' Input RDF file not found. File ',FIRDF
	stop
 198	write(*,*)' Xmol prototype file not found. File ',FXMOL
	stop
 199	write(*,*)' Torsion list not found. File ',FANGL
	stop
 921	write(*,*)' Error in renumeration array, file ',FFILTR
	stop
	end  
*          
*========================= COR3 ======================================
*
	subroutine COR3(IGR,ICL)
	include 'tranal.h'
*  calculation of 3-body correlations
	IT1=ITS(IS1(IGR))     ! num of type
	IT2=ITS(IS2(IGR))     
	IT3=ITS(IS3(IGR))    
	NBG1=ISADDR(IT1)+IS1(IGR)-ISADR(IT1)
	NEN1=ISADDR(IT1+1)
	NST1=NSITS(IT1) 
	NBG2=ISADDR(IT2)+IS2(IGR)-ISADR(IT2)
	NEN2=ISADDR(IT2+1)
	NST2=NSITS(IT2) 
	NBG3=ISADDR(IT3)+IS3(IGR)-ISADR(IT3)
	NEN3=ISADDR(IT3+1)
	NST3=NSITS(IT3) 
        NB12=(sqrt(RB12)-RMIN12)/DR12+1
	do I=NBG1,NEN1,NST1
	  do J=NBG2,NEN2,NST2
	    DX12=SX(I)-SX(J)
	    DY12=SY(I)-SY(J)
	    DZ12=SZ(I)-SZ(J)
	    call PBC(DX12,DY12,DZ12)
	    R12S=DX12**2+DY12**2+DZ12**2   
	    NR=NA*sqrt(R12S)/RDFCUT+1
	    if(IPRINT.ge.7.and.R12S.lt.RMX12)
     +write(*,*)I,J,sqrt(R12S),NR
	    if(NR.le.NA)RDF(NR)=RDF(NR)+1.
	    if(R12S.gt.RMN12.and.R12S.lt.RMX12)then
*	write(*,'(6I4,3f12.3)')I,J,IT1,IT2,IS1(IGR),IS2(IGR),RMN12,R12S,RMX12
	      R12=sqrt(R12S)
	      COR2=COR2+1.d0                   
	      do K=NBG3,NEN3,NST3
	        if(I.ne.K.and.J.ne.K)then
	          DX13=SX(I)-SX(K)
	          DY13=SY(I)-SY(K)
	          DZ13=SZ(I)-SZ(K)
	          call PBC(DX13,DY13,DZ13)
	          R13=DX13**2+DY13**2+DZ13**2
		  SCP=(DX12*DX13+DY12*DY13+DZ12*DZ13)/R12
		  X3=SCP-0.5d0*R12
		  Y3=dsqrt(R13-SCP**2)
		  IX3=(X3+RMAXX)/DXX
		  IY3=Y3/DYY            
* 	write(*,'(3I5,3f12.4,3I4)')I,J,K,R12,X3,Y3,IX3,IY3,NRX 
		  if(ICL.eq.1)then
	  	  if(IX3.le.NRX.and.IX3.gt.0.and.IY3.le.NRY)
     +            P3COR(IX3,IY3)=P3COR(IX3,IY3)+1.d0
		  else
	  	  if(IX3.le.NRX.and.IX3.gt.0.and.IY3.le.NRY)
     +            P3CO2(IX3,IY3)=P3CO2(IX3,IY3)+1.d0
		  end if 
	  	end if
	      end do
	    end if
	  end do
	end do
	return
	end
*
*============================== ORIENT ==============================
*	                                                               
	subroutine ORIENT(IGR)
	include 'tranal.h'
	dimension RV(3),D1(3),D2(3),E0(3),EX(3),EY(3),EZ(3)
	dimension FR(NS)
	data INIT/0/
	IZY=RZMAX/DDZ
*  This is only for 2-strand-DNA  (1st and 2nd types) 
	if(INIT.eq.0.and.IPROC.eq.4)then
	   INIT=1
	   NREF=0
	   XR0=0.
	   YR0=0.
	   ZR0=0.
	   do ITYP=1,2
	     IBEG=ISADR(ITYP)+1
	     IEND=ISADR(ITYP+1)
C   reference coordinates
	     do IS=IBEG,IEND
	       NREF=NREF+1
	       XR0=XR0+R(IS,1)
	       YR0=YR0+R(IS,2)
	       ZR0=ZR0+R(IS,3)
	     end do
	   end do
	   XR0=XR0/NREF
	   YR0=YR0/NREF
	   ZR0=ZR0/NREF
	   do ITYP=1,2
	     IBEG=ISADR(ITYP)+1
	     IEND=ISADR(ITYP+1)
	     do IS=IBEG,IEND
	       XLOC=R(IS,1)-XR0
	       YLOC=R(IS,2)-YR0
	       RLOC=sqrt(XLOC**2+YLOC**2)
*  define angle between -180 and 180
	       FREF=acos(XLOC/RLOC)
	       if(YLOC.lt.0.)FREF=-FREF
	       FR(IS)=FREF
	     end do 
	   end do
	   write(*,*)' ref. shift ',XR0,YR0,ZR0
	end if
*  Always:
	do IGR=1,NDUP   ! equivalent atoms "to"
	do II=1,NOI   ! equivalent atom groups "from"
	  if(IPROC.ge.5)then
	    ITL1=ITS(IO1(II))
	    ITL2=ITS(IO2(II))
	    ITL3=ITS(IO3(II))
	    do IML1=1,NSPEC(ITL1)
  	      ISHFT1=ISADDR(ITL1)+(IML1-1)*NSITS(ITL1)-ISADR(ITL1)
	      I=ISHFT1+IO1(II)
	      do IML2=1,NSPEC(ITL2)
  	        ISHFT2=ISADDR(ITL2)+(IML2-1)*NSITS(ITL2)-ISADR(ITL2)
	        I1=ISHFT2+IO2(II)
		D1(1)=SX(I1)-SX(I)
		D1(2)=SY(I1)-SY(I)
		D1(3)=SZ(I1)-SZ(I)
		call PBC(D1(1),D1(2),D1(3))
		RR=sqrt(D1(1)**2+D1(2)**2+D1(3)**2)
		if(RR.le.RB12+D12.and.RR.ge.RB12-D12)then 
		   if(IPRINT.ge.7)write(*,*)'find ',IML1,IML2,' at ',RR
		do IML3=1,NSPEC(ITL3)
  	          ISHFT3=ISADDR(ITL3)+(IML3-1)*NSITS(ITL3)-ISADR(ITL3)
	          I2=ISHFT3+IO3(II)
                  D2(1)=SX(I2)-SX(I)
                  D2(2)=SY(I2)-SY(I)
                  D2(3)=SZ(I2)-SZ(I)
	          call PBC(D2(1),D2(2),D2(3))
	  	  RR=sqrt(D2(1)**2+D2(2)**2+D2(3)**2)
		  if(RR.le.RB13+D13.and.RR.ge.RB13-D13)then 
		    XX=SX(I1)-SX(I2)
		    YY=SY(I1)-SY(I2)
		    ZZ=SZ(I1)-SZ(I2)
		    call PBC(XX,YY,ZZ)
		    RR=sqrt(XX**2+YY**2+ZZ**2)
		    if(RR.le.RB23+D23.and.RR.ge.RB23-D23)then
*  here we have selected all three atoms satisfying conditions 
*  COM of mobil coord.syst.
	  E0(1)=SX(I)
	  E0(2)=SY(I)
	  E0(3)=SZ(I)
	  if(IPRINT.ge.6)write(*,*)I,I1,I2,' ',NM(IO1(II)),
     +  NM(IO2(II)),NM(IO3(II))
 	  call VECT(D1,D2,RV) ! normal v. to 123 plane  (Z axis)
	  call NORM(RV,EZ,RN,3)
	  RV(1)=D1(1)+D2(1)                 ! mediana of 213 angle  ( X axis)
	  RV(2)=D1(2)+D2(2)
	  RV(3)=D1(3)+D2(3)
	  call NORM(RV,EX,RN2,3)            !  X axis
	  call VECT(EZ,EX,EY)               !  Y axis
	  IML=IML1
 	  if(IPROC.eq.5)call MEANSTR(II,IML,E0,EX,EY,EZ,IML2,IML3)
 	  JTYP=ITS(ISOR(IGR))
	  NBG=ISADDR(JTYP)+ISOR(IGR)-ISADR(JTYP)
	  NEN=ISADDR(JTYP+1)
	  NST=NSITS(JTYP)
	  do J=NBG,NEN,NST        ! over molecules with given site
*	    if(J.eq.ISHFT1+IO1(II).or.J.eq.ISHFT2+IO2(II).or.
*     +   ISHFT3+J.eq.IO3(II))go to 320
	    if(ISELF.eq.0.and.(NNUM(I1).eq.NNUM(J).or.NNUM(I2).eq.NNUM(J)
     + .or.NNUM(I).eq.NNUM(J)))go to 320
*	    if(J.eq.ISHFT1+IO1(II))go to 320
	    DX=SX(J)-E0(1)
	    DY=SY(J)-E0(2)
	    DZ=SZ(J)-E0(3)
	    call PBC(DX,DY,DZ)
	    RR2=DX**2+DY**2+DZ**2
            NCOR=NCOR+1
*   Calculating 3D distribution
	    if(RR2.lt.ROMAX2)then
	      RR=sqrt(RR2)
              XM=EX(1)*DX+EX(2)*DY+EX(3)*DZ       ! X-coord in molec.c.s. 
              YM=EY(1)*DX+EY(2)*DY+EY(3)*DZ       ! Y-coord in molec.c.s. 
              ZM=EZ(1)*DX+EZ(2)*DY+EZ(3)*DZ       ! Z-coord in molec.c.s.
*  calculate 3D density 
	      NXM=(XM+RXMAX)/DDX
	      NYM=(YM+RYMAX)/DDY
	      NZM=(ZM+RZMAX)/DDZ
	      if(IPROC.eq.5)then
	        if(NXM.ge.0.and.NXM.le.NOMX.and.NYM.ge.0.
     + and.NYM.le.NOMY.and.NZM.ge.0.and.NZM.le.NOMZ)then
     	         ID3(NXM,NYM,NZM)=ID3(NXM,NYM,NZM)+1
		 ICNT3=ICNT3+1
		 if(IPRINT.ge.8)write(*,*)I,J,NXM,NYM,NZM,XM,YM,ZM
		end if
	      else if(IPROC.eq.6.and.abs(abs(ZM)-ZB).le.DZB)then  
* case of Z=ZB plane            
	        NXM=(XM+RXMAX)/DDX
	        NYM=(YM+RYMAX)/DDY
		if(IPRINT.ge.7)write(*,*)I,J,NXM,NYM,NZM,XM,YM,ZM
	if(NXM.ge.0.and.NXM.le.NOMX.and.NYM.ge.0.and.NYM.le.NOMY)then
	   	  NCOR=NCOR+1
     		  POR(NXM,NYM)=POR(NXM,NYM)+1.
		end if
	      end if
	      if(IPROC.eq.7.and.abs(YM-ZB).le.DZB)then   ! Y=ZB plane         
	        NXM=(XM+RXMAX)/DDX
	        NYM=(YM+RYMAX)/DDY
		if(IPRINT.ge.7)write(*,*)I,J,NXM,NYM,NZM,XM,YM,ZM
	if(NXM.ge.0.and.NXM.le.NOMX.and.NYM.ge.0.and.NYM.le.NOMY)then
	   	  NCOR=NCOR+1
     		  POR(NXM,NYM)=POR(NXM,NYM)+1.
		end if
	      end if
	      if(IPROC.eq.8.and.abs(XM-ZB).le.DZB)then    ! X=ZB plane        
	        NXM=(XM+RXMAX)/DDX
	        NYM=(YM+RYMAX)/DDY
	if(NXM.ge.0.and.NXM.le.NOMX.and.NYM.ge.0.and.NYM.le.NOMY)then
	   	  NCOR=NCOR+1
     		  POR(NXM,NYM)=POR(NXM,NYM)+1.
		end if
	      end if    ! IPROC.eq.3
	    end if ! of (RR2.lt.
 320	    continue
	  end do    ! of J=NBG...
C----------------> ending IPROC>=5
	              end if ! RR.le.RB23
	            end if ! RR.le.RB13
	          end do  !  IML3=
		end if  ! RR.le.RB12
	      end do   ! IML2=
            end do ! IML1=
	  else      !  IPROC.ne.5
	  ITL=ITS(IO1(II))    ! type of "from" molecule
	  do IML=1,NSPEC(ITL)  ! cycle over "from" molecules
	  if(IPROC.eq.-1)then
*   This is only for DNA analysis
	     if(IO2(II).eq.0)go to 30
	     I1=IO1(II)
	     if(II.lt.10) then
 	       I2=IO1(II+1)
	       I3=IO1(II+10)
	       I4=IO1(II+11)
	     else if(II.eq.10)then
		I2=IO1(1)
		I3=IO1(20)
		I4=IO1(11)
	     else if(II.eq.11)then
		I2=IO1(20)
		I3=IO1(1)
		I4=IO1(10)
	     else
		I2=IO1(II-1)
		I3=IO1(II-10)
		I4=IO1(II-11)
	     end if
*   (0,0,0) point
	     E0(1)=0.25*(SX(I1)+SX(I2)+SX(I3)+SX(I4))
	     E0(2)=0.25*(SY(I1)+SY(I2)+SY(I3)+SY(I4))
	     E0(3)=0.25*(SZ(I1)+SZ(I2)+SZ(I3)+SZ(I4))
*   local coordinate system
             D1(1)=SX(I1)-SX(I3)
             D1(2)=SY(I1)-SY(I3)
             D1(3)=SZ(I1)-SZ(I3)
             D2(1)=SX(I2)-SX(I4)
             D2(2)=SY(I2)-SY(I4)
             D2(3)=SZ(I2)-SZ(I4)
	     if(D1(3).gt. HBOZL)D1(3)=D1(3)-BOZL
	     if(D1(3).lt.-HBOZL)D1(3)=D1(3)+BOZL
	     if(D2(3).gt. HBOZL)D2(3)=D2(3)-BOZL
	     if(D2(3).lt.-HBOZL)D2(3)=D2(3)+BOZL
	     if(IPRINT.gt.6)write(*,*)' base ',I1,I2,I3,I4
*   Other cases:
	  else if(IPROC.ge.0.and.IPROC.le.3)then
*   3 atom difining local coord.syst.   IPROC >= 0
	     ISHIFT=ISADDR(ITL)+(IML-1)*NSITS(ITL)-ISADR(ITL)
	     I=ISHIFT+IO1(II)
	     I1=ISHIFT+IO2(II)
	     I2=ISHIFT+IO3(II)
             D1(1)=SX(I1)-SX(I)
             D1(2)=SY(I1)-SY(I)
             D1(3)=SZ(I1)-SZ(I)
             D2(1)=SX(I2)-SX(I)
             D2(2)=SY(I2)-SY(I)
             D2(3)=SZ(I2)-SZ(I)
	     E0(1)=SX(I)
	     E0(2)=SY(I)
	     E0(3)=SZ(I)
	     call PBC(D1(1),D1(2),D1(3))
	     call PBC(D2(1),D2(2),D2(3))
	     if(IPRINT.ge.7)write(*,*)I,I1,I2,' ',NM(IO1(II)),
     +  NM(IO2(II)),NM(IO3(II))
C   specially for DNA - RMC calc.  DNA should be types 1,2
	  else if(IPROC.eq.4)then
C   calculating shift
	     XSA=0. 
	     YSA=0.
	     ZSA=0.
	     do ITYP=1,2
	       IBEG=ISADR(ITYP)+1
	       IEND=ISADR(ITYP+1)
	       do IS=IBEG,IEND
	         XSA=XSA+SX(IS)
		 YSA=YSA+SY(IS)
		 ZSA=ZSA+SZ(IS)
	       end do
	     end do
	     XSA=XSA/NREF
	     YSA=YSA/NREF
	     ZSA=ZSA/NREF
	     E0(1)=XSA-XR0
	     E0(2)=YSA-YR0
	     E0(3)=ZSA-ZR0
C   calculating turn around Z axis
	     FSH=0.
	     do ITYP=1,2
	       IBEG=ISADR(ITYP)+1
	       IEND=ISADR(ITYP+1)
	       do IS=IBEG,IEND
		 XREA=SX(IS)-XSA
		 YREA=SY(IS)-YSA
		 RREA=sqrt(XREA**2+YREA**2)
		 FREA=acos(XREA/RREA)
		 if(YREA.lt.0.)FREA=-FREA
		 FROT=FREA-FR(IS)
		 if(FROT.gt.PI/2.)FROT=FROT-PI
		 if(FROT.lt.-PI/2.)FROT=FROT+PI
		 FSH=FSH+FROT
	       end do
	     end do
	     FSH=FSH/NREF
	     if(IPRINT.ge.6)write(*,*)XSA,YSA,ZSA,FSH*TODGR
*  define vectors
	     EX(1)=cos(FSH)
	     EX(2)=sin(FSH)
	     EX(3)=0.
	     EY(1)=-sin(FSH)
	     EY(2)=cos(FSH)
	     EZ(3)=0.
	     EZ(1)=0.
	     EZ(2)=0.
	     EZ(3)=1.
	     go to 123
	  else
	     write(*,*)' wrong value of IPROC ',IPROC
	     stop
	  end if
	  call VECT(D1,D2,RV) ! normal v. to 123 plane  (Z axis)
	  call NORM(RV,EZ,RN,3)
	  RV(1)=D1(1)+D2(1)                 ! mediana of 213 angle  ( X axis)
	  RV(2)=D1(2)+D2(2)
	  RV(3)=D1(3)+D2(3)
	  call NORM(RV,EX,RN2,3)            !  X axis
	  call VECT(EZ,EX,EY)               !  Y axis