tranal.f

分子动力学模拟程序

 123	  if(IPROC.le.0.or.IPROC.eq.4.or.IPROC.eq.5)
     +      call MEANSTR(II,IML,E0,EX,EY,EZ,IML,IML)
 	  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.ISHIFT+IO1(II).or.J.eq.ISHIFT+IO2(II).or.
     +   ISHIFT+J.eq.IO3(II))go to 120
	    DX=SX(J)-E0(1)
	    DY=SY(J)-E0(2)
	    DZ=SZ(J)-E0(3)
*	    if(IPROC.eq.-1)then
*	      if(DZ.gt. HBOZL)DZ=DZ-BOZL 
*	      if(DZ.lt.-HBOZL)DZ=DZ+BOZL 
*	    else
	      call PBC(DX,DY,DZ)
*	    end if
	    RR2=DX**2+DY**2+DZ**2
            NCOR=NCOR+1
*   Calculating angularly average ion distribution around DNA
	    if(IPROC.eq.4)then
	       DRX=SX(J)-XSA
	       DRY=SY(J)-YSA
	       DRZ=0.
	       call PBC(DRX,DRY,DRZ)
	       RRA=sqrt(DRX**2+DRY**2)
	       NR=NA*RRA/RDFCUT+1
	       if(IPRINT.ge.7)write(*,*)I,J,RRA,NR
	       if(NR.le.NA)RDF(NR)=RDF(NR)+1.
	    end if
	    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.
	      if(IPROC.le.0.or.IPROC.eq.4.or.IPROC.eq.5)then
*  calculate 3D density 
		 NXM=(XM+RXMAX)/DDX
		 NYM=(YM+RYMAX)/DDY
		 NZM=(ZM+RZMAX)/DDZ
		 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
		   if(IPRINT.ge.7)write(*,*)I,J,NXM,NYM,NZM,XM,YM,ZM
		 end if
	      else if(IPROC.eq.1.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.2.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.3.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     !  RR2.lt.ROMAX
		  if(NCOR.le.0)then
		     write(*,*)' Integer overflow!!!'
		     stop
		  end if
 120	    continue
	    end do     !  J=NBG    
	  end do      !  IML=
 30	  continue
	  end if    ! IPROC.eq.5
	  end do
*  Mark hydrogens of the water molecule 
	  DX1=0.
	  DY1=1.
	  DZ1=0.
	  XN2=1.
	  YN2=0.
	  ZN2=0.
	  RRH2=DX1**2+DY1**2+DZ1**2
	  RRH=sqrt(RRH2)
	  CF=(DX1*XN2+DY1*YN2+DZ1*ZN2)/(RN2*RRH)
	  XM=RRH*CF                   
	  YM=sqrt(RRH2-XM**2)  
	  NXM=(XM+RXMAX)/DDX
	  NYM=YM/DDY  
	  if(IPROC.eq.1.and.XM+RXMAX.gt.0.d0.and.NXM.le.NOMX*2
     +    .and.NYM.le.NOMY)
     +    POR(NXM,NYM)=POR(NXM,NYM)+1.
	  if(IPROC.eq.2.and.NXM.ge.0.and.NXM.le.NOMX*2)
     +    POR(NXM,0)=POR(NXM,0)+1.
	  if(IPROC.eq.3.and.NYM.le.NOMY)
     +    POR(IZY,NYM)=POR(IZY,NYM)+1.
	end do
	return
	end
*
*==================== MCOR3 ==========================================
*
	subroutine MCOR3 
	include 'tranal.h'
	dimension ZI(0:NMX)
c     
	write(*,*)' Enter MCOR3'
        NB12=(sqrt(RB12)-RMIN12)/DR12+1 
	NOP3=0
	do I=1,NDUR
	  NOP3=NOP3+NSPEC(ITS(IS3(I)))
	end do
	if(LCOL2)then
	  NOP4=0
	  do I=NDUR+1,NDUP
	    NOP4=NOP4+NSPEC(ITS(IS3(I))) 
	  end do
	end if
	X0=-RMAXX             ! koord. polya vyvoda
	X1= RMAXX
	Y0=-RMAXY
	Y1= RMAXY             
	XA1=-0.5*RB12
	XA2=-XA1              ! koord. pervogo i vtorogo atomov
	YA1=0.               
	YA2=0.
	FCC=1.
*	if(IS1(1).eq.IS2(1))FCC=2.   
*   Write data file for 3-d plot (Mathematica)
	open(unit=15,file='3cor.data',status='unknown')  
 201  format(202f10.4) 
	do IYT=-NRY,NRY
	  IY=iabs(IYT)
	  Y3=IYT*DYY
	  YIY=(IY+0.5)*DYY             
	  SH3=abs(2.*PI*DXX*DYY*YIY)
	  FACR=NDUP*VOL/(COR2*NOP3*SH3*FCC) 
	  do IX=0,NRX
	    X3=IX*DXX-RMAXX     
	    ZI(IX)=P3COR(IX,IY)*FACR
	    if(IPRINT.ge.7)write(15,'(3f10.5)')X3,Y3,ZI
	  end do
	  write(15,201)(ZI(I),I=0,NRX)
	end do 
	close(15)
	write(*,*)' num. of part.',NOP3,NOP4
	do IYT=0,NTY             ! abs. coord. na grafike    
        YI=(IYT+0.5)*RMAXY/NTY 
	  IY=YI/DYY
	  YIY=(IY+0.5)*DYY             
	  SH3=abs(2.*PI*DXX*DYY*YIY)
	  FACR=NDUP*VOL/(COR2*NOP3*SH3*FCC)
	  if(LCOL2)FAC2=NDUP*VOL/(COR2*NOP4*SH3*FCC)  
	  if(IPRINT.ge.5)write(*,*)IYT,FACR,FAC2
	  do IXT=-NTX,NTX
	    XI=(IXT+0.5)*RMAXX/NTX               ! dec.coord.tochki
	    IX=(XI-X0)/DXX+0.001
	    if(IX.le.NRX.and.IX.ge.0.and.IY.le.NRY.and.IY.ge.0)then
	      COREL=P3COR(IX,IY)*FACR
	      CORE2=P3CO2(IX,IY)*FAC2
	    else
	      COREL=0.               
	      CORE2=0.
	    end if 
	    if(LCOL2)then
 	      ICOL=NCOLOR2(COREL,CORE2,AMPL,MAXCOLOR)
	      if(IPRINT.ge.6)write(*,'(2I4,2f8.3,2x,2I5,2f12.6,I5)')
     +IX,IY,XI,YI,IXT,IYT,COREL,CORE2,ICOL
	    else
   	      ICOL=NCOLOR(COREL,AMPL,MAXCOLOR) 
	      if(ICOL.gt.MAXCOLOR)ICOL=MAXCOLOR
	      if(IPRINT.ge.6)write(*,'(2I4,2f8.3,2x,2I5,f12.6,I5)')
     +IX,IY,XI,YI,IXT,IYT,COREL,ICOL
	    end if
          ITAB(IYT,IXT)=ICOL 
	  end do
	end do
*  Mark point 1,2
	IX1=XA1*NTX/RMAXX
	IX2=XA2*NTX/RMAXX
	IY =YA1*NTY/RMAXY 
	if(LCOL2)then
	  ITAB(IY,IX1)=(MAXCOLOR+1)**2
	  ITAB(IY,IX2)=(MAXCOLOR+1)**2
	  call MAKEPS2(fil3cor,IX1,IX2,IY)
	else
	  ITAB(IY,IX1)=MAXCOLOR
	  ITAB(IY,IX2)=MAXCOLOR
	  call MAKEPS(fil3cor,IX1,IX2,IY)
	end if
	return
	end
*
*=================== MORI ===============================
*
	subroutine MORI 
	include 'tranal.h'  
	real*8 TMP(0:NOMXM)
	X0=-RXMAX+DDX			! koord. polya vyvoda
	X1= RXMAX
	Y0=-RYMAX+DDY
	Y1= RYMAX
	Z0=-RZMAX+DDZ
	Z1= RZMAX
	open(unit=31,file=filori,status='unknown')
C   writing angularly averaged density distribution
	if(IPROC.eq.4)then
	  open(unit=16,file='dist.dat',status='unknown')
	  CONU=0.
	  write(16,*)'#   distribution of',NM(ISOR(1))
	  FACT=1.d27/(AVSNO*IAN)
	  do I=1,NA
	    RR=(I-0.5)*RDFCUT/NA
	    DRR=RDFCUT/NA
  	    SH12=2.*PI*DRR*RR*BOZL
	    CONU=CONU+RDF(I)/(20*IAN)
	    RDFM=RDF(I)*FACT/SH12
	    RDFU=RDF(I)*VOL/(SH12*NSPEC(ITS(ISOR(1)))*IAN)
	    if(CONU.gt.1.d-8)write(16,'(f8.3,3f12.5)')
     +      RR,RDFM,RDFU,CONU
	  end do
	  close(16)
	end if
	if(IPROC.le.0.or.IPROC.eq.4.or.IPROC.eq.5)then
*  writing 3D density in gOpenMol format
	   open(unit=31,file=filori,status='unknown')
	   N3=3
	   N200=200
	   write(31,*)N3,N200
	   write(31,*)NOMZ,NOMY,NOMX
	   write(31,'(6e13.5)')Z0,Z1,Y0,Y1,X0,X1
	   FACORI=VOL/(DDX*DDY*DDZ*NCOR)
	   write(*,*)' writing 3D distribution ...'
	   write(*,'(i8,3f8.2,3f8.4,e13.5)')
     +NCOR,RXMAX,RYMAX,RZMAX,DDX,DDY,DDZ,FACORI
	   do IZ=1,NOMZ
	      do IY=1,NOMY
		 do IX=1,NOMX
		    TMP(IX)=ID3(IX,IY,IZ)*FACORI
		 end do
		 write(31,'(6e13.5)')(TMP(IX),IX=1,NOMX)
	      end do
	   end do
	   close(31)
	      open(unit=32,file=fcrd,status='unknown')
	      write(32,'(i5,a)')NMSTR,
     &   '      average structure from MolDynamix'
	      write(32,'(a)')
     +'At        x       y        z     ...        RMS'
	      I1=1
	      do I=1,NMSTR
	        XX=XAV(I)/IAA(I)
		XX2=abs(XAV2(I)/IAA(I)-XX**2)
	        YY=YAV(I)/IAA(I)
		YY2=abs(YAV2(I)/IAA(I)-YY**2)
	        ZZ=ZAV(I)/IAA(I)
		ZZ2=abs(ZAV2(I)/IAA(I)-ZZ**2)
		RR=sqrt(XX2+YY2+ZZ2)
		write(32,'(a5,3f8.4,4x,f9.5)')NMM(I),XX,YY,ZZ,RR
	      end do
	      close(32)
	   return
	end if
	XA=0.
	YA=0. 
	NOP2=2*NOP
	NH2O=NSPEC(ITS(IO1(1)))
	NMOL=NSPEC(ITS(ISOR(1))) 
	if(DZB.ge.ZB)then 
	  SH=4.d0*DDO**2*(DZB+ZB)
	else
	  SH=8.d0*DDO**2*DZB
	end if
	FACR=NXM*NYM*1./NCOR
*	FACR=VOL/(NH2O*NMOL*SH*IAN)               
	do IY=0,NOMY
*	  do IX=0,NOMX
*	    COREL=POR(IX,IY)*FACR 
* 	    ICOL=NCOLOR(COREL,AMPL,MAXCOLOR)
*	    if(ICOL.gt.MAXCOLOR)ICOL=MAXCOLOR
*            if(IPRINT.ge.6)write(*,'(2I4,2f8.3,2x,2I5,f12.6,I5)')
*     +IX,IY,XI,YI,IXT,IYT,COREL,ICOL
*	    ITAB(IYT,IXT)=ICOL
*	  end do
	  write(31,'(500f8.4)')(POR(IX,IY)*FACR,IX=0,NOMX)
	end do
*  Mark C.O.
*	IX =XA*NTX/ROMAX
*	IY =YA*NTY/ROMAX  
*	ITAB(IX,IY)=MAXCOLOR
*	call MAKEPS(filori)
	return
	end
*
*======================== MEANSTR ====================================
* 
	subroutine MEANSTR(II,IML,E0,EX,EY,EZ,IML2,IML3)
	include 'tranal.h'
	dimension E0(3),EX(3),EY(3),EZ(3)
*  Sum over atoms which mean positions should be defined
	if(IO3(II).eq.0)return
**	do I=1,NMSTR
	I=1
 5	if(I.gt.NMSTR)go to 20
	   IIS=ICF(II)
	   ISS=ISM(I,IIS)           ! site
	   ITL=ITS(IO1(II))    ! type of "from" molecule
 	   if(IPROC.eq.5)then
 	      ITYP=ITS(ISS)
 	      if(ISREP(I).eq.-1)IML=IML2
 	      if(ISREP(I).eq.-2)IML=IML3
 	      IS=ISADDR(ITYP)+(IML-1)*NSITS(ITYP)+ISS-ISADR(ITYP)
 	   else
	     ISHIFT=ISADDR(ITL)+(IML-1)*NSITS(ITL)-ISADR(ITL)
	     IS=ISS+ISHIFT
 	   end if
	   NMM(I)=NM(NSITE(IS))(1:2)
	   DX=SX(IS)-E0(1)
	   DY=SY(IS)-E0(2)
	   DZ=SZ(IS)-E0(3)
	   if(IPROC.ne.4)call PBC(DX,DY,DZ)
	   XL=DX*EX(1)+DY*EX(2)+DZ*EX(3)
	   YL=DX*EY(1)+DY*EY(2)+DZ*EY(3)
	   ZL=DX*EZ(1)+DY*EZ(2)+DZ*EZ(3)
	   if(ISREP(I).gt.1.and.IAA(I).gt.0)then
*   check the next atom and accumulated average
	     I1=I+1
	     ISS2=ISM(I1,IIS)+ISHIFT
	     DX2=SX(ISS2)-E0(1)
	     DY2=SY(ISS2)-E0(2)
	     DZ2=SZ(ISS2)-E0(3)
	     if(IPROC.eq.0)call PBC(DX2,DY2,DZ2)
	     XL2=DX2*EX(1)+DY2*EX(2)+DZ2*EX(3)
	     YL2=DX2*EY(1)+DY2*EY(2)+DZ2*EY(3)
	     ZL2=DX2*EZ(1)+DY2*EZ(2)+DZ2*EZ(3)
	     AX0=XAV(I)/IAA(I)
	     AY0=YAV(I)/IAA(I)
	     AZ0=ZAV(I)/IAA(I)
	     AX1=XAV(I1)/IAA(I1)
	     AY1=YAV(I1)/IAA(I1)
	     AZ1=ZAV(I1)/IAA(I1)
	     RR1=(XL-AX0)**2+(YL-AY0)**2+(ZL-AZ0)**2           !  1->1
	     RR2=(XL2-AX0)**2+(YL2-AY0)**2+(ZL2-AZ0)**2        !  2->1
	     RR11=(XL-AX1)**2+(YL-AY1)**2+(ZL-AZ1)**2          !  1->2
	     RR21=(XL2-AX1)**2+(YL2-AY1)**2+(ZL2-AZ1)**2       !  2->2
	     if(ISREP(I).eq.3)then
	       I2=I+2
	       ISS3=ISM(I2,IIS)+ISHIFT
	       DX3=SX(ISS3)-E0(1)
	       DY3=SY(ISS3)-E0(2)
	       DZ3=SZ(ISS3)-E0(3)
               if(IPROC.ne.4)call PBC(DX3,DY3,DZ3)
	       AX2=XAV(I2)/IAA(I2)
	       AY2=YAV(I2)/IAA(I2)
	       AZ2=ZAV(I2)/IAA(I2)
	       XL3=DX3*EX(1)+DY3*EX(2)+DZ3*EX(3)
	       YL3=DX3*EY(1)+DY3*EY(2)+DZ3*EY(3)
	       ZL3=DX3*EZ(1)+DY3*EZ(2)+DZ3*EZ(3)
	       RR3 =(XL3-AX0)**2+(YL3-AY0)**2+(ZL3-AZ0)**2       !  3->1
	       RR31=(XL3-AX1)**2+(YL3-AY1)**2+(ZL3-AZ1)**2       !  3->2
	       RR32=(XL3-AX2)**2+(YL3-AY2)**2+(ZL3-AZ2)**2       !  3->3
	       RR12=(XL -AX2)**2+(YL -AY2)**2+(ZL -AZ2)**2       !  1->3
	       RR22=(XL2-AX2)**2+(YL2-AY2)**2+(ZL2-AZ2)**2       !  2->3
*  evaluation of different transpositions
	       ICS=1
	       R123=RR1+RR21+RR32
	       RMIN=R123
	       R213=RR2+RR11+RR32
	       if(R213.lt.RMIN)then
		  ICS=2
		  RMIN=R213
	       end if
	       R132=RR1+RR31+RR22
	       if(R132.lt.RMIN)then
		  ICS=3
		  RMIN=R132
	       end if
	       R231=RR11+RR31+RR12
	       if(R231.lt.RMIN)then
		  ICS=4
		  RMIN=R231
	       end if
	       R312=RR3+RR11+RR22
	       if(R312.lt.RMIN)then
		  ICS=5
		  RMIN=R312
	       end if
	       R321=RR3+RR21+RR12
	       if(R321.lt.RMIN)ICS=6
	       if(ICS.eq.2.or.ICS.eq.4)then       !   2 -> 1
		  XAV(I)=XAV(I)+XL2
		  YAV(I)=YAV(I)+YL2
		  ZAV(I)=ZAV(I)+ZL2
	          XAV2(I)=XAV2(I)+XL2**2
	          YAV2(I)=YAV2(I)+YL2**2
	          ZAV2(I)=ZAV2(I)+ZL2**2
	if(IPRINT.ge.7)write(*,'(3I5,a4,3f9.3)')
     +  I,IS,ICS,NM(NSITE(IS)),XL2,YL2,ZL2
	          IAA(I)=IAA(I)+1
		  if(ICS.eq.2)then              !  1->2  3
		     I1=I+1
		     I2=I+2
		  else       !  ICS=4              3->2; 1->3
		     I1=I+2
		     I2=I+1
		  end if
		  XAV(I1)=XAV(I1)+XL
		  YAV(I1)=YAV(I1)+YL
		  ZAV(I1)=ZAV(I1)+ZL
	if(IPRINT.ge.7)write(*,'(3I5,a4,3f9.3)')
     +  I,ISS2,ICS,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(I2)=XAV(I2)+XL3
		  YAV(I2)=YAV(I2)+YL3
		  ZAV(I2)=ZAV(I2)+ZL3
	if(IPRINT.ge.7)
     +write(*,'(3I5,a4,3f9.3)')I,ISS3,ICS,NM(NSITE(ISS3)),XL3,YL3,ZL3
	          XAV2(I2)=XAV2(I2)+XL3**2
	          YAV2(I2)=YAV2(I2)+YL3**2
	          ZAV2(I2)=ZAV2(I2)+ZL3**2
		  IAA(I1)=IAA(I1)+1
		  IAA(I2)=IAA(I2)+1
		  I=I+2
		  go to 10
	       else if (ICS.ge.5)then   !        3 -> 1
		  XAV(I)=XAV(I)+XL3
		  YAV(I)=YAV(I)+YL3
		  ZAV(I)=ZAV(I)+ZL3
	          XAV2(I)=XAV2(I)+XL3**2
	          YAV2(I)=YAV2(I)+YL3**2
	          ZAV2(I)=ZAV2(I)+ZL3**2
	if(IPRINT.ge.7)write(*,'(3I5,a4,3f9.3)')
     +I,IS,ICS,NM(NSITE(IS)),XL3,YL3,ZL3
	          IAA(I)=IAA(I)+1
		  if(ICS.eq.5)then                !    312
		     I1=I+1
		     I2=I+2
		  else   !   ICS=6                     321
		     I1=I+2
		     I2=I+1
		  end if
		  XAV(I1)=XAV(I1)+XL
		  YAV(I1)=YAV(I1)+YL
		  ZAV(I1)=ZAV(I1)+ZL
	if(IPRINT.ge.7)write(*,'(3I5,a4,3f9.3)')
     +I,ISS2,ICS,NM(NSITE(ISS2)),XL,YL,ZL
	          XAV2(I1)=XAV2(I1)+XL**2
	          YAV2(I1)=YAV2(I1)+YL**2
	          ZAV2(I1)=ZAV2(I1)+ZL**2