forcee.f

分子动力学模拟程序

*
      CX    = SX(L)-SX(K)
      CY    = SY(L)-SY(K)
      CZ    = SZ(L)-SZ(K)
      call PBC(AX,AY,AZ)
      call PBC(BX,BY,BZ)
      call PBC(CX,CY,CZ)
*
       AXAX = AX*AX 
       BXBX = BX*BX 
       CXCX = CX*CX
       AYAY = AY*AY 
       BYBY = BY*BY 
       CYCY = CY*CY
       AZAZ = AZ*AZ 
       BZBZ = BZ*BZ 
       CZCZ = CZ*CZ
*
       AXBX = AX*BX 
       AXCX = AX*CX 
       BXCX = BX*CX
       AXBY = AX*BY 
       AXCY = AX*CY 
       BXCY = BX*CY
       AXBZ = AX*BZ 
       AXCZ = AX*CZ 
       BXCZ = BX*CZ
*
       AYBX = AY*BX 
       AYCX = AY*CX 
       BYCX = BY*CX
       AYBY = AY*BY 
       AYCY = AY*CY 
       BYCY = BY*CY
       AYBZ = AY*BZ 
       AYCZ = AY*CZ 
       BYCZ = BY*CZ
*
       AZBX = AZ*BX 
       AZCX = AZ*CX 
       BZCX = BZ*CX
       AZBY = AZ*BY 
       AZCY = AZ*CY 
       BZCY = BZ*CY
       AZBZ = AZ*BZ 
       AZCZ = AZ*CZ 
       BZCZ = BZ*CZ
*
      DX    = AYBZ-AZBY
      DY    = AZBX-AXBZ
      DZ    = AXBY-AYBX
*
      EX    = BYCZ-BZCY
      EY    = BZCX-BXCZ
      EZ    = BXCY-BYCX
*
      D2    = DX*DX+DY*DY+DZ*DZ
      D1    = DSQRT(D2)
*
      E2    = EX*EX+EY*EY+EZ*EZ
      E1    = DSQRT(E2)
*
      DE    = DSQRT(D2*E2)
*
      DDOTE = DX*EX+DY*EY+DZ*EZ
      DII   = 1.D0/D1
      EII   = 1.D0/E1
      DEI   = DII*EII
*
      COS1F = DDOTE*DEI
      COS1F = DMIN1(COS1F, 1.D0)
      COS1F = DMAX1(COS1F,-1.D0)
      COS2F =  2.D0*COS1F*COS1F-1.D0
      COS3F = COS1F*(2.D0*COS2F-1.D0)
      COS1F2= COS1F*COS1F
      COS2F2= COS2F*COS2F
      COS3F2= COS3F*COS3F
      ARG   = DACOS(COS1F)
*
      FX1   = BZ*EY-BY*EZ
      FY1   = BX*EZ-BZ*EX
      FZ1   = BY*EX-BX*EY
*
      FX4   = BZ*DY-BY*DZ
      FY4   = BX*DZ-BZ*DX
      FZ4   = BY*DX-BX*DY
*
      FX3   =-2.D0*BX*(AYCY+AZCZ)+AX*(BYCY+BZCZ)+CX*(AYBY+AZBZ)
      FY3   =-2.D0*BY*(AXCX+AZCZ)+AY*(BXCX+BZCZ)+CY*(AXBX+AZBZ)
      FZ3   =-2.D0*BZ*(AXCX+AYCY)+AZ*(BXCX+BYCY)+CZ*(AXBX+AYBY)
*
      FX2   =-FX1-FX3 
      FY2   =-FY1-FY3 
      FZ2   =-FZ1-FZ3 
*
      GX1   = 2.D0*FX4
      GY1   = 2.D0*FY4
      GZ1   = 2.D0*FZ4
*
      GX3   = 2.D0*(AZ*DY-DZ*AY)
      GY3   = 2.D0*(AX*DZ-DX*AZ)
      GZ3   = 2.D0*(AY*DX-DY*AX)
*
      GX2   = -GX1-GX3
      GY2   = -GY1-GY3
      GZ2   = -GZ1-GZ3
*
      HX2   = 2.D0*(CZ*EY-CZ*EY)
      HY2   = 2.D0*(CX*EZ-CX*EZ)
      HZ2   = 2.D0*(CY*EX-CY*EX)
*
      HX4   = 2.D0*FX1
      HY4   = 2.D0*FY1
      HZ4   = 2.D0*FZ1
*
      UTORS = FN1 * (1.D0 + COS1F)
     X       +FN2 * (1.D0 - COS2F)
     X       +FN3 * (1.D0 + COS3F)
*
      TR(M)     = TR(M)+ARG
      TE(M)     = TE(M)+UTORS
      PINT	= PINT+UTORS
*
      FOR   = FN1 - 4.D0*FN2*COS1F + 3.D0*FN3*(2.D0*COS2F+1.D0)
      FOR   = -DEI*FOR
      FORK  = 0.5D0*COS1F*E1*DII
      FORL  = 0.5D0*COS1F*D1*EII
*
      F1X   = FOR*(FX1-FORK*GX1)
      F1Y   = FOR*(FY1-FORK*GY1)
      F1Z   = FOR*(FZ1-FORK*GZ1)
*
      F2X   = FOR*(FX2-FORK*GX2-FORL*HX2)
      F2Y   = FOR*(FY2-FORK*GY2-FORL*HY2)      
      F2Z   = FOR*(FZ2-FORK*GZ2-FORL*HZ2)
*
      F4X   = FOR*(FX4-FORL*HX4)
      F4Y   = FOR*(FY4-FORL*HY4)
      F4Z   = FOR*(FZ4-FORL*HZ4)
*
      HX(I) = HX(I)+F1X
      HY(I) = HY(I)+F1Y
      HZ(I) = HZ(I)+F1Z
*
      HX(J) = HX(J)+F2X
      HY(J) = HY(J)+F2Y
      HZ(J) = HZ(J)+F2Z
*
      HX(K) = HX(K)-F1X-F2X-F4X
      HY(K) = HY(K)-F1Y-F2Y-F4Y
      HZ(K) = HZ(K)-F1Z-F2Z-F4Z
*
      HX(L) = HX(L)+F4X
      HY(L) = HY(L)+F4Y
      HZ(L) = HZ(L)+F4Z
*
	  VIRFX  = VIRFX-(AX+BX)*FFI1-BX*FFJ1+CX*FFL1
	  VIRFY  = VIRFY-(AY+BY)*FFI2-BY*FFJ2+CY*FFL2
	  VIRFZ  = VIRFZ-(AZ+BZ)*FFI3-BZ*FFJ3+CZ*FFL3
      END DO! OF N
      END DO! OF M
*
      DO M      = NTBEG,NTEND
      TT(M)     = TT(M)+TR(M)*FNSPI
      ET(M)     = ET(M)+TE(M)*FNSPI
      END DO
*
*
*
      RETURN
      END
*
*=============== STRBND ===========================================
*
      SUBROUTINE STRBND(ITYP)
*      
* This routine calculates intramolecular forces for 
* triatomic molecules (like water)
*
      include "mdee.h"
*
      DIMENSION BS(NB),BE(NB)
*
*.... ANHARMONIC FLEXIBLE WATER MODEL:
*
      PARAMETER (AAP  = 2809.56D0 , BBP = 2809.56D0)
      PARAMETER (CCP  =  687.41D0 , DDP =   2.566D0)
      PARAMETER (EEP  = -884.63D0 , FFP = -884.63D0)
      PARAMETER (GGP  =  467.31D0 ,BBMA=1.3,BBSA=1.4  )
*
      IF(NSITS(ITYP).NE.3) STOP '!!! NOT A WATER MOLECULE'
      IF( IPOT(ITYP).GT.2) STOP '!!! IPOT.GT.2  - NOT SUPPORTED! '
*
      NBBEG    = IADB(ITYP)
      NBEND    = NBBEG+NRB(ITYP)-1
*
      I1       = NBBEG
      I2       = I1+1
      I3       = I2+1
*
      IBEG     = IADDR(ITYP)+1
      IEND     = IADDR(ITYP+1)
      NSP      = NSPEC(ITYP)
      NSS      = NSITS(ITYP)
      FNSPI    = 1.D0/DFLOAT(NSP)
      FACT     = 1.D3/AVSNO/UNITE
      FACTOR   = FACT
*
      ROH1     = RB(I1)
      ROH2     = RB(I2)
      RH1H2    = RB(I3)
*
      AAA      = AAP*FACTOR
      BBB      = BBP*FACTOR
      DDD      = DDP
      IPT      = IPOT(ITYP)
      IF(IPT.EQ.2) THEN
        AAA      = AAA/DDD**2
        BBB      = BBB/DDD**2
      END IF
*
      CCC      = CCP*FACTOR
      EEE      = EEP*FACTOR
      FFF      = FFP*FACTOR
      GGG      = GGP*FACTOR
      BBM	= BBMA
      BBS	= BBSA
*
      EFAC     =  FKCALJ*1.D+3/AVSNO/UNITE
*
      DO M     = NBBEG,NBEND
      BS(M)    = 0.D0
      BE(M)    = 0.D0
      END DO
*
      AADD     = AAA*DDD
      BBDD     = BBB*DDD
*
      ISHF     = ISADDR(ITYP)
      IBB      = ISHF+1
      JBB      = ISHF+2
      KBB      = ISHF+3
*
      IF(IPOT(ITYP).EQ.2) THEN
*
      DO N     = 1,NSP
      I        = (N-1)*NSS+IBB
      J        = (N-1)*NSS+JBB
      K        = (N-1)*NSS+KBB 
*
      AX       = SX(J)-SX(I)
      AY       = SY(J)-SY(I)
      AZ       = SZ(J)-SZ(I)
*
      BX       = SX(K)-SX(I)
      BY       = SY(K)-SY(I)
      BZ       = SZ(K)-SZ(I)
*
      CX       = SX(K)-SX(J)
      CY       = SY(K)-SY(J)
      CZ       = SZ(K)-SZ(J)
*
      call PBC(AX,AY,AZ)
      call PBC(BX,BY,BZ)
      call PBC(CX,CY,CZ)
      RRA      = AX*AX+AY*AY+AZ*AZ
      RRB      = BX*BX+BY*BY+BZ*BZ
      RRC      = CX*CX+CY*CY+CZ*CZ
      ASS      = DSQRT(RRA)
      BSS      = DSQRT(RRB)
      if(ASS.gt.BBS)write(*,*)' !!! bond length ',ASS,' for ',I,J
      if(BSS.gt.BBS)write(*,*)' !!! bond length ',BSS,' for ',I,K
      CSS      = DSQRT(RRC)
      COSA     = (AX*BX+AY*BY+AZ*BZ)/ASS/BSS
      ARG      = DACOS(COSA)
*
      SA       = ASS-ROH1
      SB       = BSS-ROH2
      SH       = CSS-RH1H2
*
      EXPA     = DEXP(-DDD*SA)
      EXPAM    = 1.D0-EXPA
      EXPB     = DEXP(-DDD*SB)
      EXPBM    = 1.D0-EXPB
*
      WB1       = AAA*EXPAM**2+EEE*SA*SH
      WB2       = BBB*EXPBM**2+FFF*SB*SH
      WB3       = CCC*SH*SH   +GGG*SA*SB
      WB        = WB1+WB2+WB3
*
      BS(I1)    = BS(I1)+ASS
      BS(I2)    = BS(I2)+BSS
      BS(I3)    = BS(I3)+CSS
*!    BS(I3)    = BS(I3)+ARG
*
      BE(I1)    = BE(I1)+WB1
      BE(I2)    = BE(I2)+WB2
      BE(I3)    = BE(I3)+WB3 
	EWB	= WB1+WB2+WB3
      PINT	= PINT+EWB
*	if(IPRINT.ge.7)write(*,'(I5,f13.5,3f10.5)')
*     +N,0.001*EWB*ENERF,ASS,BSS,CSS
*
      if(ASS.gt.BBM)then
        FB3A	= -2.*SA*AADD*DDD
      else
        FB3A	= -2.*AADD*EXPA*EXPAM
      end if
      FB3A     = FB3A   - EEE*SH-GGG*SB
      if(BSS.gt.BBM)then
        FB3B	= -2.*SB*BBDD*DDD
      else
        FB3B	= -2.*BBDD*EXPB*EXPBM
      end if
      FB3B     = FB3B   - FFF*SH-GGG*SA
      FB3C     = -2.0*CCC*SH            - EEE*SA-FFF*SB
*
      FB3A     = FB3A/ASS
      FB3B     = FB3B/BSS
      FB3C     = FB3C/CSS
*
      AXX       = FB3A*AX
      AYY       = FB3A*AY
      AZZ       = FB3A*AZ
*
      BXX       = FB3B*BX
      BYY       = FB3B*BY
      BZZ       = FB3B*BZ
*
      CXX       = FB3C*CX
      CYY       = FB3C*CY
      CZZ       = FB3C*CZ
*
      HX(I)    = HX(I)-AXX-BXX
      HY(I)    = HY(I)-AYY-BYY
      HZ(I)    = HZ(I)-AZZ-BZZ
*
      HX(J)    = HX(J)+AXX   -CXX
      HY(J)    = HY(J)+AYY   -CYY
      HZ(J)    = HZ(J)+AZZ   -CZZ
*
      HX(K)    = HX(K)   +BXX+CXX
      HY(K)    = HY(K)   +BYY+CYY
      HZ(K)    = HZ(K)   +BZZ+CZZ
*  7.3.6  Virial calculations
	VIX	= AX*AXX + BX*BXX + CX*CXX
	VIY	= AY*AYY + BY*BYY + CY*CYY
	VIZ	= AZ*AZZ + BZ*BZZ + CZ*CZZ
	VIRB	= VIRB + VIX+VIY+VIZ
	VIRFX	= VIRFX+VIX
	VIRFY	= VIRFY+VIY
	VIRFZ	= VIRFZ+VIZ
*
      END DO! OF N
*
      ELSE
*
      DO N     = 1,NSP
      I        = (N-1)*NSS+IBB
      J        = (N-1)*NSS+JBB
      K        = (N-1)*NSS+KBB
*
      AX       = SX(J)-SX(I)
      AY       = SY(J)-SY(I)
      AZ       = SZ(J)-SZ(I)
*
      BX       = SX(K)-SX(I)
      BY       = SY(K)-SY(I)
      BZ       = SZ(K)-SZ(I)
*
      CX       = SX(K)-SX(J)
      CY       = SY(K)-SY(J)
      CZ       = SZ(K)-SZ(J)
      call PBC(AX,AY,AZ)
      call PBC(BX,BY,BZ)
      call PBC(CX,CY,CZ)
*
      RRA      = AX*AX+AY*AY+AZ*AZ
      RRB      = BX*BX+BY*BY+BZ*BZ
      RRC      = CX*CX+CY*CY+CZ*CZ
*
      ASS      = DSQRT(RRA)
      BSS      = DSQRT(RRB)
      CSS      = DSQRT(RRC)
      COSA     = (AX*BX+AY*BY+AZ*BZ)/ASS/BSS
      ARG      = DACOS(COSA)
*
      SA       = ASS-ROH1
      SB       = BSS-ROH2
      SH       = CSS-RH1H2
*
      WB1      = AAA*SA*SA+EEE*SA*SH
      WB2      = BBB*SB*SB+FFF*SB*SH 
      WB3      = CCC*SH*SH+GGG*SA*SB
      WB       = WB1+WB2+WB3
*
      BS(I1)    = BS(I1)+ASS
      BS(I2)    = BS(I2)+BSS
      BS(I3)    = BS(I3)+CSS
*!    BS(I3)    = BS(I3)+ARG
*
      BE(I1)    = BE(I1)+WB1
      BE(I2)    = BE(I2)+WB2
      BE(I3)    = BE(I3)+WB3
      PINT	= PINT+WB1+WB2+WB3
*
      FB3A     = -2.0*AAA*SA - EEE*SH-GGG*SB
      FB3B     = -2.0*BBB*SB - FFF*SH-GGG*SA
      FB3C     = -2.0*CCC*SH - EEE*SA-FFF*SB
*
      FB3A     = FB3A/ASS
      FB3B     = FB3B/BSS
      FB3C     = FB3C/CSS
*
      AXX       = FB3A*AX
      AYY       = FB3A*AY
      AZZ       = FB3A*AZ
*
      BXX       = FB3B*BX
      BYY       = FB3B*BY
      BZZ       = FB3B*BZ
*
      CXX       = FB3C*CX
      CYY       = FB3C*CY
      CZZ       = FB3C*CZ
*
      HX(I)    = HX(I)-AXX-BXX
      HY(I)    = HY(I)-AYY-BYY
      HZ(I)    = HZ(I)-AZZ-BZZ
*
      HX(J)    = HX(J)+AXX   -CXX
      HY(J)    = HY(J)+AYY   -CYY
      HZ(J)    = HZ(J)+AZZ   -CZZ
*
      HX(K)    = HX(K)   +BXX+CXX
      HY(K)    = HY(K)   +BYY+CYY
      HZ(K)    = HZ(K)   +BZZ+CZZ
*
*  7.4.6 Virial calculations                      
	VIX	= AX*AXX + BX*BXX + CX*CXX
	VIY	= AY*AYY + BY*BYY + CY*CYY
	VIZ	= AZ*AZZ + BZ*BZZ + CZ*CZZ
	VIRB	= VIRB + VIX+VIY+VIZ
	VIRFX	= VIRFX+VIX
	VIRFY	= VIRFY+VIY
	VIRFZ	= VIRFZ+VIZ
      END DO! OF N
*
      END IF
*
      DO M      = NBBEG,NBEND
      BB(M)     = BB(M)+BS(M)*FNSPI
      EB(M)     = EB(M)+BE(M)*FNSPI
      END DO! OF M
*      print *,bs(1)*fnspi
*     X       ,bs(2)*fnspi
*     X       ,bs(3)*fnspi
*
	RETURN
      END
*
*=============== SLOW FORCES =============================================
*
      SUBROUTINE FORCES
*
*  LJ and Re-part of electrostatic interactions
*
      include "mdee.h"
*
      timel0	= cputime(0.d0)
      DO N         = 1,NSTOT
        FX(N)        = 0.D0