forcee.f

分子动力学模拟程序

*  Collection of force calculation
*
*=============== INTERNAL FORCES =============================================
*
*   Covalent bonds
*
*=============== GETBND ===========================================
*
      SUBROUTINE GETBND(ITYP)
*
      include "mdee.h"
* 
      character*3 NAMN
      DIMENSION BR(NB),BE(NB)
*                    
      NRBS      = NRB(ITYP)
      IF(NRBS.LE.0) RETURN
      timeb0	= cputime(0.d0)
*
      ISBEG     = ISADDR(ITYP)+1
      ISEND     = ISADDR(ITYP +1)
*     
      NAMN=NAME(ITYP)(1:3)
      NSS       = NSITS (ITYP)
      IF(IPOT(ITYP).eq.1.or.IPOT(ITYP).eq.2)THEN
        CALL STRBND(ITYP)
        timeb	= timeb+cputime(timeb0)
        RETURN
      END IF
*
      NBBEG     = IADB(ITYP)
      NBEND     = NBBEG+NRBS-1
*
      DO M      = NBBEG,NBEND
        BR(M)     = 0.D0
        BE(M)     = 0.D0
      END DO
      NSP       = NSPEC (ITYP)
      FNSPI     = 1.D0/DFLOAT(NSP)
*
      ISHF      = ISADDR(ITYP)
      DO M      = NBBEG,NBEND
        IBB       = IB(M)+ISHF
        JBB       = JB(M)+ISHF
        REQ       = RB(M)
        FCN       = FB(M)
*
        DO K      = 1,NSP
          I         = (K-1)*NSS+IBB
          J         = (K-1)*NSS+JBB
          BX        = SX(J)-SX(I)
          BY        = SY(J)-SY(I)
          BZ        = SZ(J)-SZ(I)
          call PBC(BX,BY,BZ)
          BSQ       = BX*BX+BY*BY+BZ*BZ
          BD        = DSQRT(BSQ)
*     if(ityp.eq.1.and.m.eq.nbbeg) print *,bd
          R1I       = 1.D0/BD
          BBB       =  BD-REQ
          DBND      = BBB*FCN
          EBDD	    = DBND*BBB
          FORB      =  -2.0D0*DBND*R1I
          BR(M)     = BR(M)+BD
          BE(M)     = BE(M)+EBDD
          PINT	    = PINT+EBDD
*          write(*,*)I,J,IS,JS,NM(IS),NM(JS),BD,EBDD*0.001*ENERF
          HX(I)     = HX(I)-BX*FORB
          HY(I)     = HY(I)-BY*FORB
          HZ(I)     = HZ(I)-BZ*FORB
          HX(J)     = HX(J)+BX*FORB
          HY(J)     = HY(J)+BY*FORB
          HZ(J)     = HZ(J)+BZ*FORB
	  VIRFX	    = VIRFX+FORB*BX**2  
	  VIRFY	    = VIRFY+FORB*BY**2  
	  VIRFZ	    = VIRFZ+FORB*BZ**2
        END DO! OF K
*
      END DO! OF M
*
      DO M      = NBBEG,NBEND
        BB(M)     = BB(M)+BR(M)*FNSPI
        EB(M)     = EB(M)+BE(M)*FNSPI
      END DO
*
      timeb	= timeb+cputime(timeb0)
      RETURN
      END
*
*=============== GETANG ===========================================
*
      SUBROUTINE GETANG(ITYP)
*
      include "mdee.h"
*
      DIMENSION AR(NA),AE(NA)
*
      NRAS      = NRA(ITYP)
      IF(NRAS.LE.0) RETURN
*
      timea0	= cputime(0.d0)
      NABEG     = IADA(ITYP)
      NAEND     = NABEG+NRAS-1
*
      DO M      = NABEG,NAEND
        AR(M)     = 0.D0
        AE(M)     = 0.D0
      END DO
*
      ISHF      = ISADDR(ITYP)
      IBEG      = IADDR (ITYP)+1
      IEND      = IADDR (ITYP+1)
      NSP       = NSPEC (ITYP)
      NSS       = NSITS (ITYP)
      FNSPI     = 1.D0/DFLOAT(NSP)
*
      DO M      = NABEG,NAEND
        II        = IA(M)+ISHF
        JJ        = JA(M)+ISHF
        KK        = KA(M)+ISHF
        FCN       = FA(M)
        AEQ       = RA(M)
        DO L      = 1,NSP
          I         = (L-1)*NSS+II
          J         = (L-1)*NSS+JJ
          K         = (L-1)*NSS+KK
*
          AX        = SX(I)-SX(J)
          AY        = SY(I)-SY(J)
          AZ        = SZ(I)-SZ(J)
          call PBC(AX,AY,AZ)
          ASQ       = AX*AX+AY*AY+AZ*AZ
          ASQI      = 1.D0/ASQ
*
          BX        = SX(K)-SX(J)
          BY        = SY(K)-SY(J)
          BZ        = SZ(K)-SZ(J)
          call PBC(BX,BY,BZ)
          BSQ       = BX*BX+BY*BY+BZ*BZ
          BSQI      = 1.D0/BSQ
*
          AB        = DSQRT(ASQ*BSQ)
          ABI       = 1.D0/AB
*
          COSA      = AX*BX+AY*BY+AZ*BZ
          COSA      = COSA*ABI
          COSA      = DMIN1(COSA, 1.D0)
          COSA      = DMAX1(COSA,-1.D0)
          ALF       = DACOS(COSA)
          AVA       = AVA+ALF
          SINA      = DSIN(ALF)
          SINAI     =  1.D0/SINA
          DALF      =  ALF-AEQ
          ABC       = DALF*FCN
          EANG	    = ABC*DALF	
*
          AR(M)     = AR(M)+ALF
          AE(M)     = AE(M)+EANG
          PINT	    = PINT+EANG 
*
*      if(nstep.eq.2)
*     x print  "(2i5,5x,3(1x,a4),f8.2,3x,f18.8)"
*     X ,m,l,nam(i),nam(j),nam(k)
*     x ,alf*180./3.1412,ABC*DALF*enerf*1.d-3
          SAB       =-2.D0*ABC*SINAI
          CAB       = SAB*COSA
          FAB       = SAB*ABI
          FAA       = CAB*ASQI
          FBB       = CAB*BSQI
*
          FAX       = FAB*BX-FAA*AX
          FAY       = FAB*BY-FAA*AY
          FAZ       = FAB*BZ-FAA*AZ
*
          FBX       = FAB*AX-FBB*BX
          FBY       = FAB*AY-FBB*BY
          FBZ       = FAB*AZ-FBB*BZ
*
          HX(I)     = HX(I)-FAX
          HY(I)     = HY(I)-FAY
          HZ(I)     = HZ(I)-FAZ
          HX(J)     = HX(J)+FAX+FBX
          HY(J)     = HY(J)+FAY+FBY
          HZ(J)     = HZ(J)+FAZ+FBZ
          HX(K)     = HX(K)-FBX
          HY(K)     = HY(K)-FBY
          HZ(K)     = HZ(K)-FBZ
	    VIRFX      = VIRFX-AX*FAX-BX*FBX
	    VIRFY      = VIRFY-AY*FAY-BY*FBY
	    VIRFZ      = VIRFZ-AZ*FAZ-BZ*FBZ
        END DO! OF NN
      END DO! OF M
*
      DO M      = NABEG,NAEND
      AA(M)     = AA(M)+AR(M)*FNSPI
      EA(M)     = EA(M)+AE(M)*FNSPI
      END DO
*
      timea	= timea+cputime(timea0)
      RETURN
      END
*
*=============== GETTRS ============================================== 
*
      SUBROUTINE GETTRS(ITYP)
*
      include "mdee.h"
*
      DIMENSION TR(NT),TE(NT)
*
      NRTS      = NRT(ITYP)
      IF(NRTS.LE.0) RETURN
      timet0	= cputime(0.d0)
      IF(ITORS(ITYP).NE.0) THEN
      CALL GETTOR(ITYP)
      timet	= timet+cputime(timet0)
      RETURN
      END IF
*
      NTBEG     = IADT(ITYP)
      NTEND     = NTBEG+NRTS-1
*
      DO M      = NTBEG,NTEND
        TR(M)     = 0.D0
        TE(M)     = 0.D0
      END DO
*
      ISHF      = ISADDR(ITYP)
      IBEG      = IADDR (ITYP)+1
      IEND      = IADDR (ITYP+1)
      NSP       = NSPEC (ITYP)
      NSS       = NSITS (ITYP)
      FNSPI     = 1.D0/DFLOAT(NSP)
*
      DO M      = NTBEG,NTEND
        II        = IT(M)+ISHF
        JJ        = JT(M)+ISHF
        KK        = KT(M)+ISHF
        LL        = LT(M)+ISHF
        FCN       = FT(M)
        TEQ       = RT(M)
        MUL       = NMUL(M)
        DO N      = 1,NSP
          I         = (N-1)*NSS+II
          J         = (N-1)*NSS+JJ
          K         = (N-1)*NSS+KK
          L         = (N-1)*NSS+LL
*      if(n.eq.1) print *,'torsio  ',nam(i),nam(j),nam(k),nam(l)
*
          ax        = sx(j)-sx(i)
          ay        = sy(j)-sy(i)
          az        = sz(j)-sz(i)
          bx        = sx(k)-sx(j)
          by        = sy(k)-sy(j)
          bz        = sz(k)-sz(j)
          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)
*
          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)
            si        = dsin(arg)
            if( dabs(si).lt.1.0d-12 ) si = dsign(1.0d-12,si)
*
*  For different types of torsion angles:
*
*   ---------------------------------------------
*   6.4.1  Normal AMBER-type torsion angle 
          if(ITORS(M).eq.0)then
            FCN       = FT(M)
            TEQ       = RT(M)
            MUL       = NMUL(M)
            EARG      = MUL*ARG-TEQ
            POTT      = FCN*(1.D0+DCOS(EARG))  
            DERI      = -FCN*MUL*DSIN(EARG)
*   6.4.2 MM3 force field torsional angle
	  else if(ITORS(M).eq.1)then
            FN1   = FT1(M)
            FN2   = FT2(M)
            FN3   = FT3(M)  
	    COS1F = dcos(ARG)
            COS2F =  2.D0*COS1F**2-1.D0
            COS3F = COS1F*(2.D0*COS2F-1.D0)
	    SIN2F = 2.d0*SIN1F*COS1F
	    SIN3F = SIN1F*COS2F+SIN2F*COS1F
            POTT  = FN1 * (1.D0 + COS1F)
     X             +FN2 * (1.D0 - COS2F)
     X             +FN3 * (1.D0 + COS3F)  
	    DERI	= -FN1*SIN1F+2.d0*FN2*SIN2F-3.d0*FN3*SIN3F
*   6.4.3            torsional angle 
	  else if(ITORS(M).eq.5)then  
	    FN0   = FT(M)
            FN1   = FT1(M)
            FN2   = FT2(M)
            FN3   = FT3(M)  
            FN4   = FT4(M)
            FN5   = FT5(M)
            TEQ   = RT(M)
            EARG  = ARG-TEQ
            COSI      =  DCOS(EARG)
            COSI2     =  COSI*COSI
            COSI3     =  COSI*COSI2
            COSI4     =  COSI*COSI3
            COSI5     =  COSI*COSI4
	    SINI	=  dsin(EARG)
            POTT =    FN0+ FN1*COSI+FN2*COSI2+     FN3*COSI3+
     +                              FN4*COSI4+     FN5*COSI5 
	    DERI =  -SINI*(FN1+2.d0*FN2*COSI +3.d0*FN3*COSI2+
     +                         4.d0*FN4*COSI3+5.d0*FN5*COSI4)
*   6.4.4 Improper dihedral angle
	  else if(ITORS(M).eq.-1)then                        
            FCN       = FT(M)
            TEQ       = RT(M)
            EARG      = ARG-TEQ
            POTT      = FCN*EARG**2  
            DERI      = 2.d0*FCN*EARG
	  else
	    write(*,*)'  Torsion type ',ITORS(M),' not supported'
	    write(*,*)'  Torsion ',I,J,K,L,' on ',NAME(ITYP)
	    stop
          end if
*
*...Energies:
*
            TR(M)     = TR(M)+ARG
            TE(M)     = TE(M)+POTT
            PINT	= PINT+POTT
	  if(IPRINT.ge.9)write(*,*)I,J,K,L,ARG*TODGR,POTT/EFACT
*
*...Forces:
*
            HELP      = -FCN*MUL*DSIN(EARG)
            de1       = 1.0d0/den/si*HELP
            axb       = axb/den*co
            bxc       = bxc/den*co
*X
      dnum      =  cx*bt - bx*bc
      dden      = (ab*bx - ax*bt)*bxc
      FFI1      = (dnum  - dden) * de1
      dnum      = ((bx-ax)*bc - ab*cx ) + (2.0*ac*bx - cx*bt)
      dden      = axb*(bc*cx-bx*ct) + (ax*bt-at*bx-ab*(bx-ax))*bxc
      FFJ1      = (dnum - dden) * de1
      dnum      = ab*bx - ax*bt
      dden      = axb*( bt*cx - bc*bx )
      FFL1      = (dnum - dden) * de1
      FFK1      = -(ffi1+ffj1+ffl1)
*
      HX(I)     = HX(I)+ffi1
      HX(J)     = HX(J)+ffj1
      HX(K)     = HX(K)+ffk1
      HX(L)     = HX(L)+ffl1
*Y
      dnum      =  cy*bt - by*bc
      dden      = (ab*by - ay*bt)*bxc
      FFI2      = (dnum  - dden) * de1
      dnum      = ((by-ay)*bc - ab*cy ) + (2.0*ac*by - cy*bt)
      dden      = axb*(bc*cy-by*ct) + (ay*bt-at*by-ab*(by-ay))*bxc
      FFJ2      = (dnum - dden) * de1
      dnum      = ab*by - ay*bt
      dden      = axb*( bt*cy - bc*by )
      FFL2      = (dnum - dden) * de1
      FFK2      = -(ffi2+ffj2+ffl2)
      HY(I)     = HY(I)+ffi2
      HY(J)     = HY(J)+ffj2
      HY(K)     = HY(K)+ffk2
      HY(L)     = HY(L)+ffl2
*Z
      dnum      =  cz*bt - bz*bc
      dden      = (ab*bz - az*bt)*bxc
      FFI3      = (dnum  - dden) * de1
      dnum      = ((bz-az)*bc - ab*cz ) + (2.0*ac*bz - cz*bt)
      dden      = axb*(bc*cz-bz*ct) + (az*bt-at*bz-ab*(bz-az))*bxc
      FFJ3      = (dnum - dden) * de1
      dnum      = ab*bz - az*bt
      dden      = axb*( bt*cz - bc*bz )
      FFL3      = (dnum - dden) * de1
      FFK3      = -(ffi3+ffj3+ffl3)
      HZ(I)     = HZ(I)+ffi3
      HZ(J)     = HZ(J)+ffj3
      HZ(K)     = HZ(K)+ffk3
      HZ(L)     = HZ(L)+ffl3
*   6.5.4 Contributions to projections of virial
	  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
*
          endif ! den>0
        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
*
*
      timet	= timet+cputime(timet0)
      RETURN
      END
*
*=============== GETTOR ============================================
*
      SUBROUTINE GETTOR(ITYP)
*
      include "mdee.h"
*
      DIMENSION TR(NT),TE(NT)
*
      NRTS  = NRT(ITYP)
      IF(NRTS.LE.0) RETURN
*
      NTBEG = IADT(ITYP)
      NTEND = NTBEG+NRTS-1
*
      DO M  = NTBEG,NTEND
      TR(M) = 0.D0
      TE(M) = 0.D0
      END DO
*
      ISHF  = ISADDR(ITYP)
      IBEG  = IADDR (ITYP)+1
      IEND  = IADDR (ITYP+1)
      NSP   =  NSPEC (ITYP)
      NSS   = NSITS (ITYP)
      FNSPI = 1.D0/DFLOAT(NSP)
*
      DO M  = NTBEG,NTEND
      II    = IT(M)+ISHF
      JJ    = JT(M)+ISHF
      KK    = KT(M)+ISHF
      LL    = LT(M)+ISHF
      FN1   = 0.5D0*FT1(M)
      FN2   = 0.5D0*FT2(M)
      FN3   = 0.5D0*FT3(M)
      DO N  = 1,NSP
      I     = (N-1)*NSS+II
      J     = (N-1)*NSS+JJ
      K     = (N-1)*NSS+KK
      L     = (N-1)*NSS+LL
*
      AX    = SX(J)-SX(I)
      AY    = SY(J)-SY(I)
      AZ    = SZ(J)-SZ(I)
*
      BX    = SX(K)-SX(J)
      BY    = SY(K)-SY(J)
      BZ    = SZ(K)-SZ(J)