forcee.f

分子动力学模拟程序

        FY(N)        = 0.D0
        FZ(N)        = 0.D0
        OX(N)        = 0.D0
        OY(N)        = 0.D0
        OZ(N)        = 0.D0
      END DO! OF N 
*
      IF(NTYPES.EQ.1.AND.NSPEC(1).EQ.1) RETURN
*
*
* Sum over atom pairs
      DO INP      = 1,MBLN
 10	ISP0=NBL1(INP)
	JSP=NBL2(INP)
        ISP=iabs(ISP0)
	ITYP=ITYPE(ISP)
	JTYP=ITYPE(JSP)
	ISB =NSITE(ISP)
  	JSB =NSITE(JSP)
        IMOL=NNUM(ISP)
        JMOL=NNUM(JSP)
	MTR=MDX(ITYP,JTYP)
	if(ISP0.le.0)then  
* signals that this pair is 1-4 neigbours
	  if(ITYP.ne.JTYP)write(*,*)
     +'!!! different mol. types for 1-4 interaction:',
     +ITYP,JTYP,ISP,JSP,' FF'
	  FCTLJ=C14LJ(ITYP)
*   8.2.1 special 1-4 interactions
c   Look through the list of special 1-4 interactions
C   and check whether this pair belong to it
C   This procedure may seem cumbersobe, but I didnot find better way
C   without using (NSITES,NSITES) array. Really this is not very
C   time-consuming; saving memory is preferable
	  do I=1,NNAD 
	    ILA=ILJ(I)
	    if(ILA.eq.ISB)then  
	      do J=1,NNAD 
	        JLA=ILJ(J)
	        if(JLA.eq.JSB)then
                  EPSI       =  EPAD(I)*EPAD(J)
                  SIGM       =  SIGAD(I)+SIGAD(J)
                  go to 195
	        end if
	      end do
              EPSI       =  EPAD(I)*EPSIL(JSB)
              SIGM       =  SIGAD(I)+SIGMA(JSB)
              go to 195
	    end if
	  end do
	  do J=1,NNAD
	    JLA=ILJ(J)
	    if(JLA.eq.JSB)then
              EPSI       =  EPSIL(ISB)*EPAD(J)
              SIGM       =  SIGMA(ISB)+SIGAD(J)
              go to 195
	    end if
	  end do
*   8.2.2  Normal LJ for 1-4 connected pairs
C
C   Note: EPSIL here, above and below is really sqrt(epsil) in LJ potential
C  (see also setup.f )
C
          EPSI       =  EPSIL(ISB)*EPSIL(JSB)
          SIGM       =  SIGMA(ISB)+SIGMA(JSB)
 195      continue 
          A6         = -EPSI*SIGM**6*FCTLJ    
          B12       =  EPSI*SIGM**12*FCTLJ    
  	  QIJ=CHARGE(ISB)*CHARGE(JSB)*COULF*C14EL(ITYP)
*   8.2.3 Not 1-4 neigbours
	else  ! if(ISP0
*  the same molecule - not scale in EE
          if(IMOL.eq.JMOL)then
            EPSI       =  EPS(ISB,JSB)
            SIGM       =  SIG(ISB,JSB)
            A6         = -EPSI*SIGM**6
            B12        =  EPSI*SIGM**12
            QIJ        =  CHARGE(ISB)*CHARGE(JSB)*COULF
*  scaled parameters defined in RECINT
          else
  	    A6=SIX(ISB,JSB)
	    B12=TWL(ISB,JSB)
	    QIJ=ONE(ISB,JSB)
          end if
        end if   ! if(ISPO...
        DX           = SX(ISP)-SX(JSP)
        DY           = SY(ISP)-SY(JSP)
        DZ           = SZ(ISP)-SZ(JSP)
*                  call PBC(DX,DY,DZ)
        if(DX.gt. HBOXL)DX=DX-BOXL
        if(DX.lt.-HBOXL)DX=DX+BOXL
        if(DY.gt. HBOYL)DY=DY-BOYL
        if(DY.lt.-HBOYL)DY=DY+BOYL
        if(DZ.gt. HBOZL)DZ=DZ-BOZL
        if(DZ.lt.-HBOZL)DZ=DZ+BOZL
*  
	RR           = DX**2+DY**2+DZ**2
        R1           = DSQRT(RR)
C        if(R1.le.3..and.(ISP.eq.12.or.JSP.eq.12))
C     +write(*,*)'F- long ',R1,' atoms ',ISP,JSP
        R2I          = 1.0D0/RR
	R1I	   = R1*R2I
*  Electrostatic interaction
        ALPHAR       = ALPHAD*R1
	TTT          = 1.D0/(1.D0+B1*ALPHAR)
        EXP2A        = DEXP(-ALPHAR**2)
        ERFC = ((((A5*TTT+A4)*TTT+A3)*TTT+A2)*TTT+A1)*TTT*EXP2A 
	EE1	= QIJ*R1I
*  LJ interaction
        R6I          = R2I*R2I*R2I
        ESR          = (     A6+      B12*R6I)*R6I
        FSR          = (6.D0*A6+12.D0*B12*R6I)*R6I*R2I
	if(IMOL.eq.JMOL)then 
C   This can not happen for different molecules... 
C   Scaling constant from the input 
          if(ISP0.le.0)then
	    EES=EE1
            FES=EES*R2I
          else
            FES          = QIJ*(TOTPI*ALPHAR*EXP2A+ERFC)*R1I*R2I
            EES          = EE1*ERFC
          end if
        else   ! normal case
          FES          = QIJ*(TOTPI*ALPHAR*EXP2A+ERFC)*R1I*R2I
          EES          = EE1*ERFC
	end if
        FORCE        = FES+FSR
*  Forces
          FX(ISP)      = FX(ISP)+FORCE*DX
          FY(ISP)      = FY(ISP)+FORCE*DY
          FZ(ISP)      = FZ(ISP)+FORCE*DZ
          FX(JSP)      = FX(JSP)-FORCE*DX
          FY(JSP)      = FY(JSP)-FORCE*DY
          FZ(JSP)      = FZ(JSP)-FORCE*DZ
          VIR1	 = VIR1+FSR*RR
	  if(IMOL.ne.JMOL.or.LMOVE(ITYP))then
	    VIRX	= VIRX+FORCE*DX**2  
	    VIRY	= VIRY+FORCE*DY**2  
	    VIRZ	= VIRZ+FORCE*DZ**2  
	  end if
	  POTES(MTR)   = POTES(MTR)+EE1
	  POTLJ(MTR)   = POTLJ(MTR)+ESR
	  PELS1=PELS1+EES
          PE1 = PE1+ESR
*
      END DO! OF INP					
*
*     print *,'potes(1)',potes(1)*enerf*1.d-3/256.
*     print *,'potlj(1)',potlj(1)*enerf*1.d-3/256.
*     print *,'potes(2)',potes(2)*enerf*1.d-3/256.
*     print *,'potlj(2)',potlj(2)*enerf*1.d-3/256.
*     print *,'potes(3)',potes(3)*enerf*1.d-3/256.
*     print *,'potlj(3)',potlj(3)*enerf*1.d-3/256.
*
*Accumulate total forces
*
      DO N         = 1,NSTOT
	M	= NNUM(N)
        GX(N)        = GX(N)+FX(N)
        GY(N)        = GY(N)+FY(N)
        GZ(N)        = GZ(N)+FZ(N)
	WIRS	= WIRS-GX(N)*(SX(N)-X(M))-GY(N)*(SY(N)-Y(M))-
     +               GZ(N)*(SZ(N)-Z(M))
      END DO! OF N
      timel	= timel+cputime(timel0)
*
      RETURN
      END
*
*======================= LOCALF ==========================================
*
      SUBROUTINE LOCALF
*
*  LJ and electrostatic interactions of closest neigbours
*
      include "mdee.h"
*
      times0	= cputime(0.d0)
      DO N         = 1,NSTOT
        FX(N)        = 0.D0
        FY(N)        = 0.D0
        FZ(N)        = 0.D0
      END DO! OF N
*
      IF(NTYPES.EQ.1.AND.NSPEC(1).EQ.1) RETURN
*
* Sum over atom pairs
      EEE=0.
      ELJ=0.
      DO INP      = 1,MBSH
	ISP0=NBS1(INP)
	JSP=NBS2(INP)
        ISP=iabs(ISP0)
	ITYP=ITYPE(ISP)
	JTYP=ITYPE(JSP)
	ISB =NSITE(ISP)
  	JSB =NSITE(JSP)
        IMOL=NNUM(ISP)
        JMOL=NNUM(JSP)
	MTR=MDX(ITYP,JTYP)
	if(ISP0.le.0)then  
* signals that this pair is 1-4 neigbours
	  if(ITYP.ne.JTYP)write(*,*)
     +'!!! different mol. types for 1-4 interaction:',
     +ITYP,JTYP,ISP,JSP,' FF'
	  FCTLJ=C14LJ(ITYP)
*   9.2.1 special 1-4 interactions
c   Look through the list of special 1-4 interactions
C   and check whether this pair belong to it
C   This procedure may seem cumbersobe, but I didnot find better way
C   without using (NSITES,NSITES) array. Really this is not very
C   time-consuming; saving memory is preferable
	  do I=1,NNAD 
	    ILA=ILJ(I)
	    if(ILA.eq.ISB)then  
	      do J=1,NNAD 
	        JLA=ILJ(J)
	        if(JLA.eq.JSB)then
                  EPSI       =  EPAD(I)*EPAD(J)
                  SIGM       =  SIGAD(I)+SIGAD(J)
                  go to 195
	        end if
	      end do
              EPSI       =  EPAD(I)*EPSIL(JSB)
              SIGM       =  SIGAD(I)+SIGMA(JSB)
              go to 195
	    end if
	  end do
	  do J=1,NNAD
	    JLA=ILJ(J)
	    if(JLA.eq.JSB)then
              EPSI       =  EPSIL(ISB)*EPAD(J)
              SIGM       =  SIGMA(ISB)+SIGAD(J)
              go to 195
	    end if
	  end do
*   8.2.2  Normal LJ for 1-4 connected pairs
C
C   Note: EPSIL here, above and below is really sqrt(epsil) in LJ potential
C  (see also setup.f )
C
          EPSI       =  EPSIL(ISB)*EPSIL(JSB)
          SIGM       =  SIGMA(ISB)+SIGMA(JSB)
 195      continue 
          A6         = -EPSI*SIGM**6*FCTLJ    
          B12       =  EPSI*SIGM**12*FCTLJ    
  	  QIJ=CHARGE(ISB)*CHARGE(JSB)*COULF*C14EL(ITYP)
*   9.2.3 Not 1-4 neigbours
	else  ! if(ISP0
*  the same molecule - not scale in EE
          if(IMOL.eq.JMOL)then
            EPSI       =  EPS(ISB,JSB)
            SIGM       =  SIG(ISB,JSB)
            A6         = -EPSI*SIGM**6
            B12        =  EPSI*SIGM**12
            QIJ        =  CHARGE(ISB)*CHARGE(JSB)*COULF
*  scaled parameters defined in RECINT
          else
  	    A6=SIX(ISB,JSB)
	    B12=TWL(ISB,JSB)
	    QIJ=ONE(ISB,JSB)
          end if
        end if   ! if(ISPO...
*
*   9.2 Calculate interaction parameters for the given atom pair 
*   ------------------------------------------------------------
        DX           = SX(ISP)-SX(JSP)
        DY           = SY(ISP)-SY(JSP)
        DZ           = SZ(ISP)-SZ(JSP)
*                  call PBC(DX,DY,DZ)
        if(DX.gt. HBOXL)DX=DX-BOXL
        if(DX.lt.-HBOXL)DX=DX+BOXL
        if(DY.gt. HBOYL)DY=DY-BOYL
        if(DY.lt.-HBOYL)DY=DY+BOYL
        if(DZ.gt. HBOZL)DZ=DZ-BOZL
        if(DZ.lt.-HBOZL)DZ=DZ+BOZL
*  
	RR           = DX**2+DY**2+DZ**2
        R1           = DSQRT(RR)
C        if(R1.le.3..and.(ISP.eq.12.or.JSP.eq.12))
C     +write(*,*)'F- short ',R1,' atoms ',ISP,JSP
          R2I          = 1.0D0/RR
	  R1I	   = R1*R2I
*  Electrostatic interaction
          ALPHAR       = ALPHAD*R1
	  TTT          = 1.D0/(1.D0+B1*ALPHAR)
          EXP2A        = DEXP(-ALPHAR**2)
          ERFC = ((((A5*TTT+A4)*TTT+A3)*TTT+A2)*TTT+A1)*TTT*EXP2A 
	  EE1	= QIJ*R1I
	  if(ISP0.le.0)then 
C   This can not happen for different molecules... 
            EES=EE1
            FORCE=EES*R2I
          else
            FORCE        = QIJ*(TOTPI*ALPHAR*EXP2A+ERFC)*R1I*R2I
            EES          = EE1*ERFC
	  end if
	  POTES(MTR)    = POTES(MTR)+EE1
	  PELS2=PELS2+EES
*  LJ interaction
	  if(A6.ne.0.d0.or.B12.ne.0.d0)then
            R6I          = R2I*R2I*R2I
            ESR          = (     A6+      B12*R6I)*R6I
            FSR          = (6.D0*A6+12.D0*B12*R6I)*R6I*R2I
	      FORCE        = FORCE+FSR
              VIR2		= VIR2+FSR*RR
	      POTLJ(MTR)    = POTLJ(MTR)+ESR
	      PE2		= PE2+ESR
	     end if     
*        write(*,'(4i5,5e11.4)')ISP0,JSP,ISB,JSB,R1,ESR,EES,PE
*  Forces
        FX(ISP)      = FX(ISP)+FORCE*DX
        FY(ISP)      = FY(ISP)+FORCE*DY
        FZ(ISP)      = FZ(ISP)+FORCE*DZ
        FX(JSP)      = FX(JSP)-FORCE*DX
        FY(JSP)      = FY(JSP)-FORCE*DY
        FZ(JSP)      = FZ(JSP)-FORCE*DZ
	VIRFX	= VIRFX+FORCE*DX**2  
	VIRFY	= VIRFY+FORCE*DY**2  
	VIRFZ	= VIRFZ+FORCE*DZ**2  
*
      END DO! OF INP
      PERMOL=0.001*ENERF / NOP
*
*Accumulate total forces
*
      DO N         = 1,NSTOT
	M	= NNUM(N)
        HX(N)        = HX(N)+FX(N)
        HY(N)        = HY(N)+FY(N)
        HZ(N)        = HZ(N)+FZ(N)
        WIRSS	= WIRSS-HX(N)*(SX(N)-X(M))-HY(N)*(SY(N)-Y(M))-
     +               HZ(N)*(SZ(N)-Z(M))
      END DO! OF N
      VIRS=VIRS+VIR1
      VIRLS=VIRLS+VIR1
      times	= times+cputime(times0)
*
      RETURN
      END 
*
*=============== FURIR =======================================
*
      SUBROUTINE FURIR
*
      include "mdee.h"
      real*8 TSIN(NTOT),TCOS(NTOT)
*                                       
	timef0	= cputime(0.d0)
*
      call INIFUR
*
      IF(NTYPES.EQ.1.AND.NSPEC(1).EQ.1) RETURN
	if(NKV.eq.0)return
*
      QSS         = 0.D0
      DO ITYP     = 1,NTYPES
        QSS         = QSS+QSUM(ITYP)
      END DO! OF ITYP
      IF(QSS.EQ.0.D0) RETURN
*
      QPE         = 0.D0
      CFUR=4.d0*PI*COULF/VOL     
*
      DO I        = 1,NSTOT
	  FX(I)=0.d0
	  FY(I)=0.d0
	  FZ(I)=0.d0
      END DO! OF I
*   reciprocal space Evald
      do IKV=1,NKV
        RXV=RX(IKV)
        RYV=RY(IKV)
        RZV=RZ(IKV)
        SCS=0.
        SSN=0.
C$DIR NO_RECURRENCE
        do I=1,NSTOT
          SCP=RXV*SX(I)+RYV*SY(I)+RZV*SZ(I)
          ITYP=ITYPE(I)
          IC=NSITE(I)
	  if(IEE(ITYP).eq.0)then
	      QO=CHARGE(IC)*EC(ME)
	  else
	      QO=CHARGE(IC)
	  end if
*  sin and cos from tables
          SINSC=QO*sin(SCP)
          COSSC=QO*cos(SCP)
          TSIN(I)=SINSC
          TCOS(I)=COSSC
          SSN=SSN+SINSC
          SCS=SCS+COSSC
        end do
        AKC=CFUR*AK(IKV)
        QPE=QPE+AKC*(SSN**2+SCS**2)
        AK2=2.d0*AKC
        do I=1,NSTOT
          QFOR=AK2*(SCS*TSIN(I)-SSN*TCOS(I))
          FX(I)=FX(I)+RXV*QFOR
          FY(I)=FY(I)+RYV*QFOR
          FZ(I)=FZ(I)+RZV*QFOR
        end do
	VIRX=VIRX-RXV**2*AKV*(SSN**2+SCS**2) 
	VIRY=VIRY-RYV**2*AKV*(SSN**2+SCS**2)
	VIRZ=VIRZ-RZV**2*AKV*(SSN**2+SCS**2)
      end do      
      DO I     = 1,NSTOT
        GX(I)    = GX(I)+FX(I)
        GY(I)    = GY(I)+FY(I)
        GZ(I)    = GZ(I)+FZ(I)
        IMOL	= NNUM(I)
        WIR1	= WIR1-FX(I)*(SX(I)-X(IMOL))-FY(I)*(SY(I)-Y(IMOL))-
     +               FZ(I)*(SZ(I)-Z(IMOL))
      END DO! OF I
*
      VIRX=VIRX+QPE
      VIRY=VIRY+QPE
      VIRZ=VIRZ+QPE
      PELS1=PELS1+QPE
      FACTR=  0.001*ENERF/FNOP
*	if(IPRINT.ge.7)write(*,'(a,4f12.5)')
*     +' FUR_TRUE:',FACTR*QPE
      timef	= timef+cputime(timef0)
      return
      end
*
*================= INIFUR =========================================
*
      subroutine INIFUR
	include "mdee.h"
      data INIT/0/
      save INIT
      if(INIT.ne.0)return 
        INIT=1
	if(ALPHAD.eq.0.d0)then
	  NKV=0
	  return
      end if    
      ALP2=ALPHAD**2
      CUTK2=4.d0*ALP2*FEXP
	CUTK=sqrt(CUTK2)
      FKX=2.d0*PI/BOXL
      FKY=2.d0*PI/BOYL
      FKZ=2.d0*PI/BOZL
      IKV=0
	KMAX=CUTK/FKX+1
	KMAY=CUTK/FKY+1
	KMAZ=CUTK/FKZ+1
      do IKZ=-KMAZ,KMAZ
        do IKX=0,KMAX
          if(IKX.eq.0)then
            if(IKZ.gt.0)then
              KY1=0
            else
              KY1=1
            end if
            KY2=KMAY
          else
            KY1=-KMAY
            KY2=KMAY
          end if
          do IKY=KY1,KY2
            RK2=(IKX*FKX)**2+(IKY*FKY)**2+(IKZ*FKZ)**2
            if(RK2.le.CUTK2)then
      	      IKV=IKV+1
              if(IKV.gt.NKVMAX)stop 'increase NKVMAX'
              RX(IKV)=IKX*FKX
              RY(IKV)=IKY*FKY
              RZ(IKV)=IKZ*FKZ
              AK(IKV)=exp(-0.25d0*RK2/ALP2)/RK2
            end if