forcee.f

分子动力学模拟程序

          end do
        end do
      end do      
      NKV=IKV
      write(*,*)'Ewald sum:  Kmax=',KMAX,'   Num. of k-vectors ',NKV
      return	
      end 
*   12.  Calculation of self-interaction term in the Ewald sum
*   ----------------------------------------------------------
C
C   Really this procedure calculates two terms:
C   1) Standard self-interaction Ewald term, giving only constant 
C   contribution to the energy
C   2) "Molecular correction" which is due to the fact, that bound
C   atoms in a molecule do not interact electrostatically, but 
C   reciprocal Ewald term includes these interactions also.
*
*============= ETERMS =================================================
*
      SUBROUTINE ETERMS
*  12.1 Front matter
      include "mdee.h"         
      real*8 SPEEI(NTPS)
      data INIT/0/  
      save SPEEI,SPEI,ABC,INIT
      timee0	= cputime(0.d0)
      if(LMNBL)return
      WIR1=0.
*
*  12.2  Electrostatic self-interaction
*  ------------------------------------
	INIT=1
        ABC = ALPHAD/DSQRT(PI)
	SPEI=0.
        do ITYP=1,NTYPES 
	  SPET=0.
	  do I=1,NSITS(ITYP)
	    IS=ISADR(ITYP)+I
  	    QII=ONE(IS,IS)
            SPET         = SPET-NSPEC(ITYP)*ABC*QII
	  end do         
	  SPET=SPET/NUMTASK
	  SPEEI(ITYP)=SPET
	  SPEI=SPEI+SPET
	end do
      do I=1,NTYPES
	SELFPE(I)=SPEEI(I)
      end do            
      SPE=SPEI
*
*  12.3  Calculate intramolecular correction
*  -----------------------------------------
*  12.3.1 Start cycle over bound atoms 
C  list of bound atoms is prepeared in BNDLST (setup.f)   
      IBEG=NUMTASK-TASKID 
      do ISP=IBEG,NSTOT,NUMTASK
	NSP = NNBB(ISP)
	do J=1,NSP
	  JSP=    INBB(ISP,J)
	  if(JSP.le.0)JSP=-JSP
	  if(ISP.gt.JSP)then  ! avoid double count!
	    IS	= NSITE(ISP)
	    JS	= NSITE(JSP)
	    ITYP  = ITYPE(ISP)
  	    QIJ = ONE(IS,JS)
            DX          = SX(ISP)-SX(JSP)
            DY          = SY(ISP)-SY(JSP)
            DZ          = SZ(ISP)-SZ(JSP)
*  12.3.2  Apply periodic boundary conditions
            if(DX.gt. HBOXL)DX=DX-BOXL
            if(DX.lt.-HBOXL)DX=DX+BOXL
	    if(LHEX)then
C  hexagonal periodic cell)
	      XY=BOXYC*DX
            if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then
              DY=DY-BOYL
              DX=DX-HBOXL
	      XY=BOXYC*DX
            end if
          if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then
              DY=DY+BOYL
              DX=DX+HBOXL
	      XY=BOXYC*DX
            end if
            if(DY.gt.BOXY3+XY)then
              DY=DY-BOYL
              DX=DX+HBOXL
	      XY=BOXYC*DX
            end if
            if(DY.lt.-BOXY3+XY)then
              DY=DY+BOYL
              DX=DX-HBOXL
            end if
	    else                   
C  rectangular cell
              if(DY.gt. HBOYL)DY=DY-BOYL
              if(DY.lt.-HBOYL)DY=DY+BOYL
	    end if
            if(DZ.gt. HBOZL)DZ=DZ-BOZL
            if(DZ.lt.-HBOZL)DZ=DZ+BOZL
	    if(LOCT)then
	      CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL)
	      DX=DX-sign(CORSS,DX)
	      DY=DY-sign(CORSS,DY)
	      DZ=DZ-sign(CORSS,DZ)
	    end if
*	  call PBC(DX,DY,DZ)
*   12.3.3 Calculate contributions to energy and forces
            R2          = DX**2+DY**2+DZ**2
            R2I         = 1.D0/R2
            R1          = DSQRT(R2)
            R1I         = R1*R2I
            ALPHAR      = ALPHAD*R1
*  	          call RERFC(ALPHAR,ERFC,EXP2A)
            TTT         = 1.D0/(1.D0+B1*ALPHAR)
            EXP2A       = DEXP(-ALPHAR**2)
            ERFC = ((((A5*TTT+A4)*TTT+A3)*TTT+A2)*TTT+A1)*TTT*EXP2A
            EES  = QIJ*(ERFC-1.d0)*R1I
            SPE         = SPE+EES
	    SELFPE(ITYP)  = SELFPE(ITYP)+EES
              FES  = QIJ*((TOTPI*ALPHAR*EXP2A+ERFC)-1.d0)*R1I*R2I
              GX(ISP)     = GX(ISP)+FES*DX
              GY(ISP)     = GY(ISP)+FES*DY
              GZ(ISP)     = GZ(ISP)+FES*DZ
              GX(JSP)     = GX(JSP)-FES*DX
              GY(JSP)     = GY(JSP)-FES*DY
              GZ(JSP)     = GZ(JSP)-FES*DZ
	      VIRX	= VIRX+FES*DX**2  
	      VIRY	= VIRY+FES*DY**2  
	      VIRZ	= VIRZ+FES*DZ**2
	  end if
	end do 
      END DO! OF I
 100  timee	= timee+cputime(timee0)
      RETURN
      END
*                 
*============= ELRCLJ =================================================
*
      SUBROUTINE ELRCLJ
*
*   long-range corrections to LJ energy and pressure
*
      include"mdee.h"
      data INIT/0/
*
      DO I        = 1,MOLINT
      PELRC(I)    = 0.D0
      VRLRC(I)    = 0.D0
      END DO! OF I
*
      RC          = DSQRT(RSSQ)
      DO 400 ITYP = 1    ,NTYPES
        ISBEG       = ISADR(ITYP)+1
        ISEND       = ISADR(ITYP +1)
        NSPI        = NSPEC(ITYP)
	if(IEE(ITYP).eq.0)then
	  ESCI=EC(ME)
	else
	  ESCI=1.d0
	end if
        DO 400   IS = ISBEG,ISEND
          DO 300 JTYP = ITYP ,NTYPES
            MT          = MDX  (ITYP,JTYP)
            JSBEG       = ISADR(JTYP)+1
            JSEND       = ISADR(JTYP +1)
            NSPJ        = NSPEC(JTYP)
            FNOPIJ      = DFLOAT(NSPI*NSPJ)
            if(ITYP.eq.JTYP)FNOPIJ=0.5*FNOPIJ
	    if(IEE(JTYP).eq.0)then
	      ESCJ=EC(ME)
	    else
	      ESCJ=1.d0
	    end if
            ESCL=ESCI*ESCJ
            DO 300   JS = JSBEG,JSEND
	      EP          = EPS(IS,JS)
              SI          = SIG(IS,JS)  
              PELRC(MT)   = ESCL*PELRC(MT)+EP*SI**3*
     X(( 1.D0/9.D0)*(SI/RC)**9 - ( 1.D0/3.D0)*(SI/RC)**3)*FNOPIJ
              VRLRC(MT)   = ESCL*VRLRC(MT)+EP*SI**3*
     X((-4.D0/3.D0)*(SI/RC)**9 + (      2.D0)*(SI/RC)**3)*FNOPIJ
  300     CONTINUE
  400 CONTINUE
*
      FCV	= UNITP/(3.*VOL)
      PELR	= 0.d0
      VIRD	= 0.d0     
      DO I      = 1,MOLINT
        PELRC(I)  = 4.*PI*PELRC(I)/VOL
        VRLRC(I)  = -4.*PI*VRLRC(I)/VOL
        PELR=PELR+PELRC(I)
        VIRD=VIRD+VRLRC(I)    
        if(INIT.eq.0.and.IPRINT.ge.6)then
          PRINT "('*** PE  LRC ',I2,': ',F12.8,' kJ/mol')",
     +    I,PELRC(I)*PERMOL
        PRINT "('*** PR  LRC ',I2,': ',F12.4,' atm   ')",I,VRLRC(I)*FCV
        end if
      END DO
      INIT=1
*
      RETURN
      END
*
C===========================================================================
*
*  15. Recalculation of list of neighbours
*  ---------------------------------------
C  This procedure calculate lists of closest and far neighbours
C  for subroutines LOCALF and FORCES
C  These lists must be recalculated after some time
C  Each node calculate goes through its own list of atom pairs 
C  and creates own lists, which are used in subroutines for force calculation
C
C================ CHCNB =======================================================
C     
      subroutine CHCNB(LRDF)
*  15.1 Front matter     
      include "mdee.h"
      logical LGRC,LRDF 
      parameter (RMX2=5.**2,MMOL=5)
      data INIT/0/ 
      timev0	= cputime(0.d0)
      MBLN=0
      MBSH=0
      IOVER=0
      IESH=0
      LGRC=LGR.and.LRDF
      RDFCUT2=RDFCUT**2
*
*  15.2  Organize cycle over ALL atom pairs
*  ----------------------------------------
      IBEG=NUMTASK-TASKID 
C  cycle over I,J is organized is follows:
C  it takes I>J if I and J are of the same parity and I<J if I and J have
C  different parity. This allows nearly uniformly distribute atom pairs 
C  between the nodes
      do I=1,NSTOT
	IMOL=NNUM(I) 
	ITYP=ITYPE(I)
*   15.2.1 Even strips
*        do J=I-2,1,-2
        do J=1,I-1
	  JMOL=NNUM(J)
	  DX=SX(I)-SX(J)
	  DY=SY(I)-SY(J)
	  DZ=SZ(I)-SZ(J)
*   15.2.1.1 Periodic boundary conditions
C	    call PBC (DX,DY,DZ)
	  if(DX.gt. HBOXL)DX=DX-BOXL
          if(DX.lt.-HBOXL)DX=DX+BOXL
          if(DZ.gt. HBOZL)DZ=DZ-BOZL
          if(DZ.lt.-HBOZL)DZ=DZ+BOZL
	    if(LHEX)then
C  hexagonal periodic cell)
	      XY=BOXYC*DX
            if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then
              DY=DY-BOYL
              DX=DX-HBOXL
	      XY=BOXYC*DX
            end if
          if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then
              DY=DY+BOYL
              DX=DX+HBOXL
	      XY=BOXYC*DX
            end if
            if(DY.gt.BOXY3+XY)then
              DY=DY-BOYL
              DX=DX+HBOXL
	      XY=BOXYC*DX
            end if
            if(DY.lt.-BOXY3+XY)then
              DY=DY+BOYL
              DX=DX-HBOXL
            end if
	    else                   
C  rectangular cell
            if(DY.gt. HBOYL)DY=DY-BOYL
            if(DY.lt.-HBOYL)DY=DY+BOYL
	    end if
	    if(LOCT)then
C  truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL
	      CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL)
	      DX=DX-sign(CORSS,DX)
	      DY=DY-sign(CORSS,DY)
	      DZ=DZ-sign(CORSS,DZ)
	    end if                  
	    RR2=DX**2+DY**2+DZ**2  
*  15.2.1.2 Case of "molecular" cutoff - calculate COM distance
C  If we do not use Ewald summation (LMNBL=.true.), 
C  it is important that at least
C  small molecules (less than MMOL atoms) where inside or outside
C  cut-off radius as a whole 
	    if(LMNBL)then
	      NMI=NSITS(ITYP) 
	      NMJ=NSITS(ITYPE(J))
	      if(NMI.gt.MMOL.and.NMJ.gt.MMOL)then
		RRM=RR2
		go to 110
	      else if(NMI.gt.MMOL.and.NMJ.le.MMOL)then
	        DX=SX(I)-X(JMOL)
	        DY=SY(I)-Y(JMOL)
	        DZ=SZ(I)-Z(JMOL)
	      else if(NMJ.gt.MMOL)then
	        DX=X(IMOL)-SX(J)
	        DY=Y(IMOL)-SY(J)
	        DZ=Z(IMOL)-SZ(J)
	      else 
	        DX=X(IMOL)-X(JMOL)
	        DY=Y(IMOL)-Y(JMOL)
	        DZ=Z(IMOL)-Z(JMOL)
	      end if
*	      call PBC (DX,DY,DZ)
	      if(DX.gt. HBOXL)DX=DX-BOXL
              if(DX.lt.-HBOXL)DX=DX+BOXL
              if(DZ.gt. HBOZL)DZ=DZ-BOZL
              if(DZ.lt.-HBOZL)DZ=DZ+BOZL
	      if(LHEX)then
C  hexagonal periodic cell)
	        XY=BOXYC*DX
            if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then
                DY=DY-BOYL
                DX=DX-HBOXL
	        XY=BOXYC*DX
              end if
          if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then
                DY=DY+BOYL
                DX=DX+HBOXL
	        XY=BOXYC*DX
              end if
              if(DY.gt.BOXY3+XY)then
                DY=DY-BOYL
                DX=DX+HBOXL
	        XY=BOXYC*DX
              end if
              if(DY.lt.-BOXY3+XY)then
                DY=DY+BOYL
                DX=DX-HBOXL
              end if
	      else                   
C  rectangular cell
                if(DY.gt. HBOYL)DY=DY-BOYL
                if(DY.lt.-HBOYL)DY=DY+BOYL
	      end if
	      if(LOCT)then
C  truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL
	        CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL)
	        DX=DX-sign(CORSS,DX)
	        DY=DY-sign(CORSS,DY)
	        DZ=DZ-sign(CORSS,DZ)
	      end if                  
C  square distance betwen COM of small molecules
	      RRM=DX**2+DY**2+DZ**2 
	    else   ! .not.LMNBL
	      RRM=RR2              
	    end if
*  15.2.1.3 Check if the atoms are bound
 110	    if(IMOL.eq.JMOL)then    ! bound atoms should not be in the list
	      if(RR2.le.RMX2)then
C   Check that the atoms are non-bound. Check using binding list of the
C   atom which has less bound atoms, only for atoms with distance < RMX
	        if(NNBB(I).le.NNBB(J))then
		  do K=1,NNBB(I) 
C   1-3 bound:  not in the list
		    if(INBB(I,K).eq.J)go to 125
                    if(INBB(I,K).eq.-J.and.L14NB(ITYP))then
C   1-4 bound:  put into list with negative first index
C   (these atom pairs may have a special way of force calculations)
  	              MBSH=MBSH+1
	              if(MBSH.gt.NBSMAX)IESH=1
                      if(IESH.eq.0)then
	                NBS1(MBSH)=-I
		        NBS2(MBSH)=J
                      end if
*            write(*,'(a,2i5,f12.5,2i8)')'1-4s',I,J,sqrt(RR2),MBSH,MBLN
	              go to 125
		    end if  
	  	  end do
	        else
		  do K=1,NNBB(J) 
C   1-3 bound
		    if(INBB(J,K).eq.I)go to 125
                    if(INBB(J,K).eq.-I.and.L14NB(ITYP))then
C   1-4 bound                                     
	              MBSH=MBSH+1
	              if(MBSH.gt.NBSMAX)IESH=1
                      if(IESH.eq.0)then
	                NBS1(MBSH)=-I
		        NBS2(MBSH)=J
                      end if
*            write(*,'(a,2i5,f12.5,2i8)')'1-4s',I,J,sqrt(RR2),MBSH,MBLN
                      go to 125
		    end if  
	  	  end do
	        end if
	      end if ! (RR2.le.RMAX2
	    end if   ! (IMOL.eq.JMOL
*   15.2.1.4 Pick up non-bound atoms to the list of neigbours
	    if(RRM.le.SHORT2)then
	      MBSH=MBSH+1
	      if(MBSH.gt.NBSMAX)IESH=1
              if(IESH.eq.0)then
	        NBS1(MBSH)=I
	        NBS2(MBSH)=J
              end if
*            write(*,'(a,2i5,f12.5,2i8)')'no_s',I,J,sqrt(RR2),MBSH,MBLN
	    else if(RRM.le.RSSQ)then
	      MBLN=MBLN+1
	      if(MBLN.gt.NBLMAX)then
                if(IOVER.eq.0)write(*,*)
     +     '!! list of far neigbours overfilled'
                IOVER=1
              end if
              if(IOVER.eq.0)then
	        NBL1(MBLN)=I
	        NBL2(MBLN)=J
*            write(*,'(a,2i5,f12.5,2i8)')'no_l',I,J,sqrt(RR2),MBSH,MBLN
              end if
	    end if
*   15.2.1.5 RDF calculation
C  omitted            
C   Only at this point we go through all the atom pairs
C   RDF calculation
 125        continue
          IF(LGRC.and.RR2.le.RDFCUT2.and.ME.eq.1.and.NHIST.ge.IHIST)THEN	!=====> LGR
	      ISB	= NSITE(J)
	      JSB	= NSITE(I)
*              write(*,*)ISB,JSB,RR2
	      if(NGRI(ISB).le.NGRI(JSB))then
		do JS=1,NGRI(ISB) 
		  if(IGRI(ISB,JS).eq.JSB)then 
		    IP = MGRI(ISB,JS)
	            go to 137  
		  end if
		end do
	      else
		do JS=1,NGRI(JSB)
		  if(IGRI(JSB,JS).eq.ISB)then
		    IP = MGRI(JSB,JS)
	            go to 137
		  end if  
		end do         
	      end if  
	      go to 138
 137	      R1	= sqrt(RR2)
              MM        = R1*DRDFI+1
            if(MM.le.MAX)IRDF(MM,IP) = IRDF(MM,IP)+1
*            write(*,'(2i4,f10.4,4i5)')I,J,R1,MM,IP,NRDF,IRDF(MM,IP)
	    end if ! of LGRC
 138        continue
        end do
        end do
        if(ME.eq.1.and.NHIST.ge.IHIST)NRDF=NRDF+1
        if(IOVER.ne.0)then
          write(*,*)' Insrease NBLMAX in prcm.h for this job to ',MBLN
          I=MBLN/NTOT+1
          write(*,*)' (change  to:  NBLMAX =',I,'* NTOT)'
          stop
        end if
        timev	= timev+cputime(timev0)
	if(IPRINT.lt.7)return
	write(*,*)' Task ',TASKID,
     +'total num.of interactions: ',MBSH,' and ',MBLN
        if(IESH.ne.0)then
          write(*,*)'!!! list of closest neigbours overfilled ',MBSH
          stop
        end if
	return
	end