tranal.f

分子动力学模拟程序

		  XAV(I2)=XAV(I2)+XL2
		  YAV(I2)=YAV(I2)+YL2
		  ZAV(I2)=ZAV(I2)+ZL2
	if(IPRINT.ge.7)
     +  write(*,'(3I5,a4,3f9.3)')I,ISS3,ICS,NM(NSITE(ISS3)),XL2,YL2,ZL2
	          XAV2(I2)=XAV2(I2)+XL2**2
	          YAV2(I2)=YAV2(I2)+YL2**2
	          ZAV2(I2)=ZAV2(I2)+ZL2**2
		  IAA(I1)=IAA(I1)+1
		  IAA(I2)=IAA(I2)+1
		  I=I+2
		  go to 10
	       end if     ! (RR2.lt.RR1.and.RR2.lt.RR3
	     else if((RR1+RR21).gt.(RR2+RR11)) then  !   ISREP=2,  1 <-> 2
	       I1=I+1
	       XAV(I1)=XAV(I1)+XL
	       YAV(I1)=YAV(I1)+YL
	       ZAV(I1)=ZAV(I1)+ZL
	       if(IPRINT.ge.7)
     + write(*,'(2I5,a4,3f9.3)')I1,ISS2,NM(NSITE(ISS2)),XL,YL,ZL
	       XAV2(I1)=XAV2(I1)+XL**2
	       YAV2(I1)=YAV2(I1)+YL**2
	       ZAV2(I1)=ZAV2(I1)+ZL**2
	       XAV(I)=XAV(I)+XL2
	       YAV(I)=YAV(I)+YL2
	       ZAV(I)=ZAV(I)+ZL2
	if(IPRINT.ge.7)write(*,'(2I5,a4,3f9.3)')
     +I,IS,NM(NSITE(IS)),XL2,YL2,ZL2
	       XAV2(I)=XAV2(I)+XL2**2
	       YAV2(I)=YAV2(I)+YL2**2
	       ZAV2(I)=ZAV2(I)+ZL2**2
	       IAA(I)=IAA(I)+1
	       IAA(I1)=IAA(I1)+1
	       I=I+1
	       go to 10
	     end if   !  ISEP.eq.3
	   end if     !  ISEP.ne.1
	   XAV(I)=XAV(I)+XL
	   YAV(I)=YAV(I)+YL
	   ZAV(I)=ZAV(I)+ZL
	   if(IPRINT.ge.7)write(*,'(2I5,a4,3f9.3)')
     +I,IS,NM(NSITE(IS)),XL,YL,ZL
	   XAV2(I)=XAV2(I)+XL**2
	   YAV2(I)=YAV2(I)+YL**2
	   ZAV2(I)=ZAV2(I)+ZL**2
	   IAA(I)=IAA(I)+1
 10	   I=I+1
	   go to 5
*	end do
 20	   return
	end
*          
*========================= DIFFUS ======================================
*
      subroutine DIFFUS
      include 'tranal.h'
      data IREC/1/,IINIT/0/,JINIT/0/
*  calculation of diffusion coefficients 
      ITYP		= IDF
      N           = ISADDR(ITYP)
      I           = IADDR(ITYP)
      IML=NSPEC(ITYP)
	XCM=0.
	YCM=0.
	ZCM=0.
      if(IML.gt.MAXMOL)IML=MAXMOL
      DO J        = 1,IML   ! sum over molecules
	NSADR=ISADDR(ITYP)+NSITS(ITYP)*(J-1)+IAT
        I           = I+1
        XX       = SX(NSADR)
        YY       = SY(NSADR)
        ZZ       = SZ(NSADR)
**Calculate C.O.M. vectors relative to the first atom:
**        DO IS       = NSBEG,NSEND
**          N           = N+1
**          XX       = XX+MASS(IS)*SX(N)
**          YY       = YY+MASS(IS)*SY(N)
**          ZZ       = ZZ+MASS(IS)*SZ(N)
**        END DO! OF IS
        XX       = XX+ISHX(J)*BOXL
        YY       = YY+ISHY(J)*BOYL
        ZZ       = ZZ+ISHZ(J)*BOZL
*  PBC**(-1)
	  if(IINIT.ne.0.or.IREC.ne.1)then
	    IREC1=IREC-1
	    if(IREC1.le.0)IREC1=IREC1+MTR 
	    if(XX-TRCX(IREC1,J).gt. HBOXL)then
	      ISHX(J)=ISHX(J)-1
	      XX=XX-BOXL
	    end if
	    if(XX-TRCX(IREC1,J).lt.-HBOXL)then
	      ISHX(J)=ISHX(J)+1
	      XX=XX+BOXL
	    end if
	    if(YY-TRCY(IREC1,J).gt. HBOYL)then
	      ISHY(J)=ISHY(J)-1
	      YY=YY-BOYL
	    end if
	    if(YY-TRCY(IREC1,J).lt.-HBOYL)then
	      ISHY(J)=ISHY(J)+1
	      YY=YY+BOYL
	    end if
	    if(ZZ-TRCZ(IREC1,J).gt. HBOZL)then
	      ISHZ(J)=ISHZ(J)-1
	      ZZ=ZZ-BOZL
	    end if
	    if(ZZ-TRCZ(IREC1,J).lt.-HBOZL)then
	      ISHZ(J)=ISHZ(J)+1
	      ZZ=ZZ+BOZL
	    end if
	  end if	  
*  remember 
	  TRCX(IREC,J)=XX	  
	  TRCY(IREC,J)=YY	  
	  TRCZ(IREC,J)=ZZ  
	  JM=J
	  XCM=XCM+XX
	  YCM=YCM+YY
	  ZCM=ZCM+ZZ
      END DO! OF J
      XCM=XCM/IML
      YCM=YCM/IML
      ZCM=ZCM/IML
	if(JINIT.eq.0)then
	   XCM0=XCM
	   YCM0=YCM
	   ZCM0=ZCM
	   JINIT=1
	end if
      TTIM(IREC)=FULTIM*1.d-12            ! in  s
*      if(IPRINT.ge.6)write(*,'(4f10.3,I7,2x,3I3)')
*     +FULTIM,TRCX(IREC,ICNT),TRCY(IREC,ICNT),TRCZ(IREC,ICNT),
*     +IREC,ISHX(ICNT),ISHY(ICNT),ISHZ(ICNT)
	if(IPRINT.ge.6)write(*,'(f12.1,3f12.6)')
     +    FULTIM,XCM-XCM0,YCM-YCM0,ZCM-ZCM0
*  calculate <R(t)-R(0)>
	IRB=IREC  
	IRB1=IRB+1
	if(IRB1.gt.MTR)IRB1=1
	IROUND=0
	do ITM=1,NTT
	  TBACK=FULTIM*1.d-12-ITM*DTT
*	  write(*,*)ITM,FULTIM,TBACK*1.d12,IRB,TTIM(IRB)*1.d12,IRB1,
*     + TTIM(IRB1)*1.d12
 10	  if(TTIM(IRB).le.TBACK)then
	    if(TTIM(IRB1)-TTIM(IRB).lt.5.*DTT)then
	    AL=(TBACK-TTIM(IRB))/(TTIM(IRB1)-TTIM(IRB))      
	    if(AL.lt.0.d0.or.AL.gt.1.d0)write(*,'(a,f9.5,I7,3e14.5)')
     +  ' strange ALPHA: ',AL,IRB,TTIM(IRB),TBACK,TTIM(IRB1)  
          DO J        = 1,IML   ! sum over molecules
	      X0=TRCX(IRB,J)*(1.-AL)+TRCX(IRB1,J)*AL
	      Y0=TRCY(IRB,J)*(1.-AL)+TRCY(IRB1,J)*AL
	      Z0=TRCZ(IRB,J)*(1.-AL)+TRCZ(IRB1,J)*AL
	      RAX=(X0-TRCX(IREC,J))**2
	      RAY=(Y0-TRCY(IREC,J))**2
	      RAZ=(Z0-TRCZ(IREC,J))**2
	      RR=RAX+RAY+RAZ
	      RCUMX(ITM)=RCUMX(ITM)+RAX
	      RCUMY(ITM)=RCUMY(ITM)+RAY
	      RCUMZ(ITM)=RCUMZ(ITM)+RAZ
	      RCUM(ITM)=RCUM(ITM)+RR
	      NCUM(ITM)=NCUM(ITM)+1  
	      if(IPRINT.ge.8)write(*,'(2I6,2f12.6)')
     +      ITM,J,ITM*DTT*1.d12,sqrt(RR)
	    end do   
	    end if
	  else
	    IRB1=IRB
	    IRB=IRB-1
	    if(IRB.le.0)then
	      if(IINIT.eq.0)go to 30
	      IRB=MTR
	      IROUND=1
	    end if
	    if(IRB.le.IREC.and.IROUND.eq.1)go to 30
	    go to 10
	  end if
	end do ! of ITM
 30	IREC=IREC+1
	if(IREC.gt.MTR)then
	  IREC=1
	  IINIT=1
	end if
	return
	end   
*
*==================== RESTIM ==========================================
*                                                                      
	subroutine RESTIM(IMD)
	include 'tranal.h'
	real*4 TIIN(NRT1M,NRT2M),TIINI(NRT1M,NRT2M),TIINO(NRT1M,NRT2M)
	integer IRS(NRT1M,NRT2M)
	data IINIT/0/
	if(IINIT.eq.0)then
	  IINIT=1
	  do I=1,NRT1M
	    do J=1,NRT2M
	      IRS(I,J)=0
	    end do
	  end do 
	end if 
	if(IBREAK.ne.0)then
	  do I=1,NRT1
	    do J=1,NRT2
	      IRS(I,J)=0
	    end do
	  end do 
	  write(*,*)' residence time tables set to zero'
	end if     
*   over first set of sites
	do I=1,NRT1
	  ISP=IRS1(I)
	  ITYP=ITYPE(ISP)
	  IMOL=NNUM(ISP)
	  if(LIST(ITYP).eq.0)return
*   over second set of sites
	  do J=1,NRT2
	    JSP=IRS2(J) 
	    JMOL=NNUM(JSP)
	    if(IMOL.ne.JMOL)then
	    JTYP=ITYPE(JSP)
	    if(LIST(JTYP).eq.0)return
	    if(ISP.eq.JSP)go to 90
	    DX=SX(ISP)-SX(JSP)
	    DY=SY(ISP)-SY(JSP)
	    DZ=SZ(ISP)-SZ(JSP)
	    call PBC(DX,DY,DZ)
	    RR=sqrt(DX**2+DY**2+DZ**2)
*  executed at the last call
	    if(IMD.eq.1)then
*  register all pairs which are still in contact
	      if(mod(IRS(I,J),3).ne.0)then
	        TIMR=(FULTIM*1.d-12-TIINI(I,J))*1.d12
	        TIMRS2=TIMRS2+TIMR     ! register residence time 
	        NR=TIMR*NITT/TIMM+1
		write(*,*)' add ',TIMR,' ps  for ',I,J,' R=',RR
	        if(NR.gt.0.and.NR.le.NITT)DIST2(NR)=DIST2(NR)+1 
	        ITIR2=ITIR2+1
	        if(TIMR.gt.TTER*1.d12)then
	          TIMRS1=TIMRS1+TIMR
	 	  ITIR1=ITIR1+1
		end if
	        TIMRS=TIMRS+TIMR     
	        if(NR.gt.0.and.NR.le.NITT)DISTT(NR)=DISTT(NR)+1 
	        ITIR=ITIR+1
	      end if
	      go to 90
	    end if    ! IMD.eq.1
	    if(IPRINT.ge.8)write(*,*)I,J,ISP,JSP,RR
*  standard criteria (exit during short time not accounted)
 15	    if(mod(IRS(I,J),3).eq.0)then 
*  bulk (0) -> enter hydrated shell (2)
	      if(RR.gt.RRN.and.RR.lt.RRX)then
		IRS(I,J)=IRS(I,J)+2
	 	TIINI(I,J)=FULTIM*1.d-12          ! register entrance time
	        if(IPRINT.ge.6)write(*,*)' new st',I,J
	      end if   
*  Temporarly out (1)
	    else if(mod(IRS(I,J),3).eq.1)then 
*  Out more than cut-off time (TTER) -> bulk (0) 
	      if(FULTIM*1.d-12-TIINO(I,J).gt.TTER)then
	        IRS(I,J)=IRS(I,J)-1
	        if(IPRINT.ge.6)write(*,*)' out st',I,J
	        TIMR=(TIINO(I,J)-TIINI(I,J))*1.d12
	        TIMRS2=TIMRS2+TIMR     ! register residence time 
	        NR=TIMR*NITT/TIMM+1
	        if(NR.gt.0.and.NR.le.NITT)DIST2(NR)=DIST2(NR)+1 
	        ITIR2=ITIR2+1
	        if(TIMR.gt.TTER*1.d12)then
	          TIMRS1=TIMRS1+TIMR
	 	  ITIR1=ITIR1+1
		end if
	 	go to 15
	      end if
	      if(RR.gt.RRN.and.RR.lt.RRX)IRS(I,J)=IRS(I,J)+1
*  hydrated (2)
	    else
	      if(RR.le.RRN.or.RR.ge.RRX)then
*  mark as temporarly out (1); register exit time 
	        IRS(I,J)=IRS(I,J)-1
	        TIINO(I,J)=FULTIM*1.d-12
            end if
	    end if	    
*  alt. criteria (just register in - out )	     
	    if(IRS(I,J).le.2)then
* simple in
	      if(RR.gt.RRNX.and.RR.lt.RRXN)then
	        IRS(I,J)=IRS(I,J)+3
	        TIIN(I,J)=FULTIM*1.d-12 
	        if(IPRINT.ge.6)write(*,*)' new alt',I,J
	      end if
	    else
* simple out    
	      if(RR.lt.RRNN.or.RR.gt.RRXX)then 
	        TIMR=(FULTIM*1.d-12-TIIN(I,J))*1.d12
		IRS(I,J)=IRS(I,J)-3
	        TIMRS=TIMRS+TIMR     
	        NR=TIMR*NITT/TIMM+1
	        if(NR.gt.0.and.NR.le.NITT)DISTT(NR)=DISTT(NR)+1 
	        ITIR=ITIR+1
	        if(IPRINT.ge.6)write(*,*)' out alt',I,J
	      end if
	    end if
*  Count hydrated pairs
	    if(mod(IRS(I,J),3).eq.2)IHP=IHP+1
 90	    continue
	  end if ! (IMOL.ne.JMOL
	  end do
	end do
	if(IPRINT.ge.7)write(*,*)fultim,ITIR,TIMRS/ITIR
      return
	end
*
*=================== MDIFF ============================================
*                                                                      
	subroutine MDIFF(fildif)
	include 'tranal.h'
	character*64 fildif
	open(unit=9,file=fildif,status='unknown')
	write(*,*)
     +'translational time -correlation function of type',IDF
	write(9,'(a)')' D units is 10**-9 m**2/s '
	write(9,'(a,I3)')
     +'#  translational time -correlation function of type',IDF
	write(9,*)'  t(ps)  Dav    Ddif       R(A)  Dx     Dy     Dz'
	write(*,*)NTT 
	TSHO=0.
	do I=1,NTT
	  TIM=I*DTT*1.d12
	  if(NCUM(I).gt.0)then
	    TCOR=10.*RCUM(I)/(6.d0*TIM*NCUM(I))
	    TCX=10.*RCUMX(I)/(2.d0*TIM*NCUM(I))
	    TCY=10.*RCUMY(I)/(2.d0*TIM*NCUM(I))
	    TCZ=10.*RCUMZ(I)/(2.d0*TIM*NCUM(I))
	    TSH2=RCUM(I)/NCUM(I)
	    TSH=sqrt(TSH2)
	    TSD=1.d-11*TSH*(TSH-TSHO)/(3.*DTT) 
	    TSHO=TSH
            write(9,'(f10.3,2f10.5,f10.3,3f10.5,i9)')
     +TIM,TCOR,TSD,TSH,TCX,TCY,TCZ,NCUM(I)
      write(*,'(f10.3,2f10.5,f10.3,3f10.5,i9)')
     +TIM,TCOR,TSD,TSH,TCX,TCY,TCZ,NCUM(I)
	  end if
	end do
	return
	end                                               
*
*=================== TAKETOR ============================================
*                                                                        
	subroutine TAKETOR
	include "tranal.h"
*  Torsion angle distributions
	do ITYP=1,NAG     !  type of torsion
	  do N=IADT(ITYP)+1,IADT(ITYP+1)
        I         = IT(N)  
        J         = JT(N)
        K         = KT(N)
        L         = LT(N)
	  MTYP = ITS(I)                 !  type of molecule
	  do IMOL=1,NSPEC(MTYP)         !  over molecules of this type
	    IST = ISADDR(MTYP)+(IMOL-1)*NSITS(MTYP)-ISADR(MTYP)+I
	    JST = ISADDR(MTYP)+(IMOL-1)*NSITS(MTYP)-ISADR(MTYP)+J
	    KST = ISADDR(MTYP)+(IMOL-1)*NSITS(MTYP)-ISADR(MTYP)+K
	    LST = ISADDR(MTYP)+(IMOL-1)*NSITS(MTYP)-ISADR(MTYP)+L
        ax        = sx(jst)-sx(ist)
        ay        = sy(jst)-sy(ist)
        az        = sz(jst)-sz(ist)
        bx        = sx(kst)-sx(jst)
        by        = sy(kst)-sy(jst)
        bz        = sz(kst)-sz(jst)
        cx        = sx(lst)-sx(kst)
        cy        = sy(lst)-sy(kst)
        cz        = sz(lst)-sz(kst)
*
        call PBC(AX,AY,AZ)
        call PBC(BX,BY,BZ)
        call PBC(CX,CY,CZ)
*
        ab        = ax*bx + ay*by + az*bz
        bc        = bx*cx + by*cy + bz*cz
        ac        = ax*cx + ay*cy + az*cz
        at        = ax*ax + ay*ay + az*az
        bt        = bx*bx + by*by + bz*bz
        ct        = cx*cx + cy*cy + cz*cz

        axb       = (at*bt)-(ab*ab)
        bxc       = (bt*ct)-(bc*bc)
*
        fnum      = (ab*bc)-(ac*bt)
        den       = axb*bxc
        if(den.gt.0.0d0) then 
*
          den       = dsqrt(den)         
Cosine of angle:
          co        = fnum/den
          CO        = DMIN1(CO, 1.D0)
          CO        = DMAX1(CO,-1.D0)
* Sign of angle:
      signum    = ax*(by*cz-cy*bz)+ay*(bz*cx-cz*bx)+az*(bx*cy-cx*by)
* Value of angle:
          arg       = dsign(dacos(co),signum)
	  else
	    arg=0.d0
	  end if    
	  NARG=0.5*(arg+PI)*NAA/PI+1
	  if(abs(arg*180/PI).le.120)then
	    IGS(ITYP)=IGS(ITYP)+1
	  else
	    ITRAN(ITYP)=ITRAN(ITYP)+1
	  end if
	  NPOP(ITYP,NARG)=NPOP(ITYP,NARG)+1 
	  if(IPRINT.ge.8)write(*,'(6i5,f10.3)')
     +ITYP,IMOL,IST,JST,KST,LST,ARG*180/PI
	  end do
	  end do
	end do
	end
*
*=================== MAKEPS ============================================
*                                                                       
	subroutine MAKEPS(filps,IX1,IX2,IY)
	include 'tranal.h'     
      character*2 color(0:maxcolor)
	character*12 filps
	character*80 com 
	data color/
     *   '0C', !  'bC', ! esli zadat' maxcolor==20,
     *   '1C', !  'cC', ! to nado ubrat' vosklicatel'nye
     *   '2C', !  'dC', ! znaki mezhdu dannymi
     *   '3C', !  'eC',
     *   '4C', !  'fC',
     *   '5C', !  'gC',
     *   '6C', !  'hC',
     *   '7C', !  'iC',
     *   '8C', !  'jC',
     *   '9C', !  'kC',
     *   'aC' /                        
	write(*,*)' Enter MAKEPS' 
	ICHAN=20
      open( ICHAN, file='tmp.ps', status='unknown' )
      write(20,*)IY,IX1,'(',NM(IS1(1)),') l-show'
      write(20,*)IY,IX2,'(',NM(IS2(1)),') l-show'
      do j=0,NTY
        m = -NTX
        IC = ITAB(j,1)
        do i=-NTX+1,NTX
          n = ITAB(j,i)
c n - eto uzhe No. tsveta ot 0 do maxcolor
          if( n .ne. IC ) then  
            write( 20, 100 )   j,    m, i-1, color(ic)
            IC = n
            m = i
          endif
        enddo
        write( 20, 100 )   j,    m, NTY, color(ic)
        write( 20, * )
      enddo
C      write(20,'(a)')'0 se