mdstepe.f

分子动力学模拟程序

*=============== SULTAN ================================================
*
      SUBROUTINE SULTAN
*
*  Constrained MD
*
      include "mdee.h"
*
      LOGICAL LRDF
      DATA LRDF/.false./
*
*>----------------------------
      CALL GETCOM     ! supl.f      
      call CHCNB(LRDF)       ! mdstep.f
      DO ITYP = 1,NTYPES
      CALL GETBND(ITYP)   ! i-forces.f
      CALL GETANG(ITYP)   ! i-forces.f
      CALL GETTRS(ITYP)   ! i-forces.f
      END DO!  OF ITYP 
      call LOCALF         ! l-forces.f  > not needed because RSHORT=0
*

      CALL ZEROFS(2)  ! supl.f
      CALL FURIR      ! l-forces.f
      CALL ETERMS     ! l-forces.f
      CALL FORCES     ! l-forces.f
      if(LNPT)call ELRCLJ     ! l-forces.f
      LRDF=.true.
      NSTEP=0
 1    NSTEP=NSTEP+1
      MSTEP=NSTTOT+NSTEP
*>----------------------------
*  remember old positions
*
      DO I         = 1,NSTOT
        OX(I)        = SX(I)
        OY(I)        = SY(I)
        OZ(I)        = SZ(I)
        FX(I) = GX(I)+HX(I)
        FY(I) = GY(I)+HY(I)
        FZ(I) = GZ(I)+HZ(I)
*Unconstrained velocities at step n+1/2
        VX(I)        = VX(I)+TSTEP*FX(I)
        VY(I)        = VY(I)+TSTEP*FY(I)
        VZ(I)        = VZ(I)+TSTEP*FZ(I)
*Unconstrained  postions  at step n+1
        SX(I)        = SX(I)+TSTEP*VX(I)*MASSDI(I)
        SY(I)        = SY(I)+TSTEP*VY(I)*MASSDI(I)
        SZ(I)        = SZ(I)+TSTEP*VZ(I)*MASSDI(I)
      END DO! OF L
*
*  Constrain dynamics                           
      VIRA(1)=0.
      VIRA(2)=0.
      VIRA(3)=0.
      do IMOL=1,NOP
	CALL SHEJK(IMOL)
      end do
*  first time - correct velocities
      if(MSTEP.eq.1)then
        call GETEMP
        FTMP        = SQRT(TRTEMP/TEMP)
        DO     N    = 1,NSTOT
          VX(N)       = VX(N)*FTMP
          VY(N)       = VY(N)*FTMP
          VZ(N)       = VZ(N)*FTMP
        end do
      end if
*  Realization of specified ensemble
      CALL SCALING(3)
      call GETCOM
      if(mod(MSTEP,ICHNB).eq.0)call CHCNB(LRDF)

        CALL ZEROFS(1)
	call ZEROAV
	LCSRDF=.true.
*       
*
      DO ITYP = 1,NTYPES
      CALL GETBND(ITYP)   ! i-forces.f
      CALL GETANG(ITYP)   ! i-forces.f
      CALL GETTRS(ITYP)   ! i-forces.f
      END DO!  OF ITYP 
      call LOCALF         ! l-forces.f
*
      CALL ZEROFS(2)  ! supl.f
      CALL FURIR      ! l-forces.f
      CALL ETERMS     ! l-forces.f
      CALL FORCES     ! l-forces.f  - must be the last for comp. WIRS
      if(LNPT)call ELRCLJ     ! l-forces.f
*
      timeg0	= cputime(0.d0)
      CALL SCLFRC(DONE,3)
*
*  Change of ensemble
      if(mod(MSTEP,NCEE).eq.0)call CHENS
*  Collecting averages
      CALL GETKIN     ! supl.f
*      CALL GETROT     ! supl.f
*      CALL DIPOLE
      if(mod(MSTEP,IAVER).eq.0)CALL GETAVR(1)
*
      NSTEPA=NSTEPA+1
      PEST=PELS1+PELS2+SPE
      PE = PE1+PE2+PEST
C   Total intermolecular energy
        PENL=PE*PERMOL
C   Electrostatic potential energy      
        PEL=PEST*PERMOL
C   Intramolecular potential energy      
        PENS=PINT*PERMOL
        PEKI=TKE*0.001*ENERF                 !   kinetic energy
        ETOT=PENL+PENS+PEKI                  !   total energy
        POTEN	= PENL+PENS                   !  total potential energy
        PHD=POTLJ(MOLINT)*1.d-3*ENERF/NHDM      !  energy of hydrogen bonds
        WIRSUM=WIRS+WIRSS
C   These are diffrent contribution to pressure:	
        TRYCKE	= PEST*UNITP/(3.*VOL)    ! electrostatic contr.
        TRYCKL	= (VIR1+VIR2)*UNITP/(3.*VOL)     ! LJ
        TRYCKB	= VIRB*UNITP/(3.*VOL)    ! bond (flexible)
        TRYCKA	= (VIRA(1)+VIRA(2)+VIRA(3))*UNITP/(3.*VOL) ! shejk-correction 
        TRYCKS	= WIRSUM*UNITP/(3.*VOL) ! orientational correction
        TRYCKD	= VIRD*UNITP/(3.*VOL) ! long-range correction
        DENS	= TOTMAS*1.d-3/(VOL*UNITL**3) !  current density
*
      if(IPRINT.ge.5.and.mod(MSTEP,IAVER).eq.0)
     +WRITE(6,991) MSTEP,ME,POTEN,PENS,PEL,ETOT,TEMP,TRYCKM,DENS
*  990 FORMAT('%',A1,I5,2(1X,F10.3),1X,F9.2,1X,F9.3,1X,F10.3,1X,F8.4)
          if(IPRINT.ge.6)write(*,'(2f9.6,5f10.2)')
     +   SC,SCL,TRYCK,TTR,TROT,TINT,TEMPV
	  if(IPRINT.ge.6.and.LSEP)write(*,'(4f11.1,3f11.3)')
     + TRYCKX,TRYCKY,TRYCKZ,(TRYCKX+TRYCKY+TRYCKZ)/3.,BOXL,BOYL,BOZL
          if(IPRINT.ge.6.and.IEXT.eq.1)write(*,'(2(a,e14.6))')
     +    ' E abs: ',(EABSS+EABS)*PERMOL
 	  if(IPRINT.ge.7)
     +  WRITE(6,'(8f10.2)')TRYCKL,TRYCKE,TRYCKB,TRYCKA,TRYCKS,TRYCKD
  991 FORMAT(I7,I4,4F9.3,2(1X,f9.2),f9.4)
**
      IF(mod(MSTEP,NAVER).eq.0)call GETAVR(2)
      if(LDMP.and.mod(NSTTOT+NSTEP,NDUMP).eq.0)call RESTRT(2)
      timeg	= timeg+cputime(timeg0)
*
*<---------------------
      if(NSTEP.lt.NSTEPS)go to 1
*<---------------------
*
      RETURN
      END
*
*=============== DOUBLE ================================================
*
      SUBROUTINE DOUBLE
*
      include "mdee.h"
*
      LOGICAL LRDF
      DATA LRDF/.false./
*
*>----------------------------
*    Restore forces, COM and list of neigbours
*
      LRDF=.false.
      NFREQ2=NFREQ/2+1
      VIRA(1)=0.
      VIRA(2)=0.
      VIRA(3)=0.
      CALL ZEROFS(1)
      CALL GETCOM     ! supl.f 
      call CHCNB(LRDF)      ! mdstep.f
      LRDF=.true.
      NFR2	= NFREQ/2+1
**       
      LCSRDF=.false.
      DO ITYP = 1,NTYPES
        CALL GETBND(ITYP)   ! i-forces.f
        CALL GETANG(ITYP)   ! i-forces.f
        CALL GETTRS(ITYP)   ! i-forces.f
      END DO!  OF ITYP
      call LOCALF         ! l-forces.f
      CALL SCLFRC(DONE,1)
*
*
      CALL ZEROFS(2)  ! supl.f 
      CALL FURIR      ! l-forces.f
      CALL ETERMS     ! l-forces.f
      CALL FORCES     ! l-forces.f
      call ELRCLJ     ! l-forces.f
      CALL SCLFRC(DONE,2)
*
      if(IPRINT.ge.7)then
	PELL=PELS1+PELS2+SPE
        PENL=(PE1+PE2+PELL)*PERMOL
        PEL=PELL*PERMOL
        PENS=PINT*PERMOL
        PEKI=TKE*PERMOL
        ETOT=PENL+PENS+PEKI
        WRITE(6,'(4(a6,f12.5))')
     +' Eex= ',PENL,' Ein= ',PENS,' Eel= ',PEL,' Etot=',ETOT 
      end if
*
*   Begin MD
*
      NSTEP=0
 1    NSTEP=NSTEP+1
	MSTEP=NSTTOT+NSTEP
*>----------------------------
      call SCALING(1)                 
*  velocity correction due to slow forces
      DO I           =  1,NSTOT
        VX(I)          =   VX(I)+TSTEP*GX(I)*0.5D0
        VY(I)          =   VY(I)+TSTEP*GY(I)*0.5D0
        VZ(I)          =   VZ(I)+TSTEP*GZ(I)*0.5D0
      END DO! OF I
*
*- - - - - - - FAST MOTION - - - - - - - -
*     
      LCSRDF=.false.
      call ZEROAV
      DO NSTEB       =  1,NFREQ 
	if(NSTEB.eq.NFREQ)LCSRDF=.true.
        DO I           =  1,NSTOT
          VX(I) =   VX(I)+TSTEB*HX(I)*0.5D0         ! V(t+1/2)
          VY(I) =   VY(I)+TSTEB*HY(I)*0.5D0
          VZ(I) =   VZ(I)+TSTEB*HZ(I)*0.5D0
          SX(I) =   VX(I)*TSTEB*MASSDI(I)+SX(I)     ! X(t+1)
          SY(I) =   VY(I)*TSTEB*MASSDI(I)+SY(I)
          SZ(I) =   VZ(I)*TSTEB*MASSDI(I)+SZ(I)
        END DO! OF I
*  Realization of the specified ensemble
	CALL GETCOM 
	if(NSTEB.eq.NFR2.and.mod(MSTEP,ICHNB).eq.0)then
           call CHCNB(LRDF)
        end if
*  Recalculation of slow forces
	if(mod(NSTEB,NFREQ).eq.0)then
          CALL ZEROFS(2)      ! supl.f
          CALL FURIR          ! l-forces.f
          CALL ETERMS         ! l-forces.f
          CALL FORCES         ! l-forces.f
*          if(LNPT)call ELRCLJ         ! l-forces.f
          CALL SCLFRC(DONE,2)   ! supl.f
        end if
        CALL ZEROFS(1)      ! supl.f
	DO ITYP = 1,NTYPES
          CALL GETBND(ITYP)     ! i-forces.f           F(t+1)
          CALL GETANG(ITYP)     ! i-forces.f
          CALL GETTRS(ITYP)     ! i-forces.f
        END DO!  OF ITYP
	call LOCALF            ! l-forces.f
*
        CALL SCLFRC(DONE,1)    ! supl.f
*         
          DO I           =  1,NSTOT
            VX(I)          =  VX(I)+TSTEB*HX(I)*0.5D0     !  V(t+1)
            VY(I)          =  VY(I)+TSTEB*HY(I)*0.5D0
            VZ(I)          =  VZ(I)+TSTEB*HZ(I)*0.5D0
          END DO! OF I        
*
*  Realization of the specified ensemble
*  Change of subensemble
	if(NSTEB.eq.NFR2.and.mod(MSTEP,NCEE).eq.0)call CHENS
*
        if(IPRINT.ge.8)then
          CALL GETEMP     ! aver.f
  990 FORMAT('%',A1,I5,2(1X,F10.3),1X,F9.2,1X,F10.3,1X,F9.2,1X,F9.2)
          WRITE(6,990) WHICH,NSTEB,POTEN,ETOT,TEMP,VIR,TRL,TRYCK
        end if
      END DO! OF NSTEB
*
*- - - - - - - SLOW MOTION - - - - - - - -
*     
      timeg0	= cputime(0.d0)
*
      DO I           = 1,NSTOT
        VX(I)          =  VX(I)+TSTEP*GX(I)*0.5D0
        VY(I)          =  VY(I)+TSTEP*GY(I)*0.5D0
        VZ(I)          =  VZ(I)+TSTEP*GZ(I)*0.5D0
      END DO! OF I
      CALL SCALING(2)       ! mdstep.f
*
*  Collecting averages
      NSTEPA=NSTEPA+1
      CALL GETKIN       ! supl.f
*      CALL GETROT       ! supl.f
*      CALL DIPOLE       ! supl.f
      if(mod(MSTEP,IAVER).eq.0)CALL GETAVR(1)
        PEST = PELS1+PELS2+SPE
        PE = PE1+PE2+PEST
        PENL = PE*PERMOL                ! intermolecular energy
        PEL= PEST*PERMOL                ! electrostatic energy          
        PENS=PINT*PERMOL                ! intramolecular energy
        PEKI=TKE*0.001*ENERF            ! kinetic energy
        POTEN=PENL+PENS                 ! total potential energy
        ETOT =POTEN+PEKI                ! total energy   
        WIRSUM=WIRS+WIRSS
        TRYCKE	= PEL*UNITP/(3.*VOL*PERMOL)     ! electrostatic
        TRYCKL	= (VIR1+VIR2)*UNITP/(3.*VOL)    ! LJ
        TRYCKB	= VIRB*UNITP/(3.*VOL)           ! bonds
        TRYCKS	= WIRSUM*UNITP/(3.*VOL)         ! molecular
        TRYCKA	= (VIRA(1)+VIRA(2)+VIRA(3))*UNITP/(3.*VOL) ! molecular
        TRYCKD	= VIRD*UNITP/(3.*VOL)           ! long-range corr.
        DENS	= TOTMAS*1.d-3/(VOL*UNITL**3)
*
*      if(IPRINT.ge.6)WRITE(6,'(4(a6,f12.5))')
*     +' Eex= ',PENL,' Ein= ',PENS,' Eki= ',PEKI,' Etot=',ETOT 
      if(IPRINT.ge.5.and.mod(MSTEP,IAVER).eq.0)
     +WRITE(6,991) MSTEP,ME,POTEN,PENS,PEL,ETOT,TEMP,TRYCK,DENS
          if(IPRINT.ge.6)write(*,'(2f7.4,f9.2,4f7.1,I7,I9)')
     + SC,SCL,TRYCKM,TTR,TROT,TINT,TEMPV,MBSH,MBLN   
   	  if(IPRINT.ge.6.and.LSEP)write(*,'(4f11.1,3f11.3)')
     + TRYCKX,TRYCKY,TRYCKZ,(TRYCKX+TRYCKY+TRYCKZ)/3.,BOXL,BOYL,BOZL
          if(IPRINT.ge.6.and.IEXT.eq.1)write(*,'(2(a,e14.6))')
     +    ' E abs: ',(EABS+EABSS)*PERMOL,' Ext.f. ',EFIELD
 	  if(IPRINT.ge.7)
     +  WRITE(6,'(8f10.1)')TRYCKL,TRYCKE,TRYCKB,TRYCKS,TRYCKD,TRID
  991 FORMAT(I7,I4,4F9.3,2(1X,f9.2),f9.4)
*
      if(mod(MSTEP,NAVER).eq.0)CALL GETAVR(2)
      if(LDMP.and.mod(MSTEP,NDUMP).eq.0)call RESTRT(2)
      timeg	= timeg+cputime(timeg0)
*
*<---------------------
      if(NSTEP.lt.NSTEPS)go to 1
*<---------------------
*
      RETURN
      END
*
*=============== SCALING ==============================================
*
      SUBROUTINE SCALING(IAL)
*
*  3  Realisation of the specified ensemble
*  ----------------------------------------
C  Case NVT/NPT - Nose-Hoover thermostat, "modular-invariant" version 
C  of G.J.Martyna, D.J.Tobais and M.L.Klein,
C  J.Chem.Ohys., v.101(5), p.4177 (1994)
C  
C  NVE - if specified, do simple velocity scaling if temperature 
C        deviate from reference value more that TDELT 
C
C  IAL = 1 - first half-step in double time step algorithm
C  IAL = 2 - second half-step in double time step algorithm
C  IAL = 3 - Leap-frog algorithm
C
*
*  3.1 Definitions
*  ---------------
      include "mdee.h"   
      parameter (FP3=0.5/3.)
      real*8 DPE(4),DTE(NTPS)
      data DKEO/0./,ISKK/0/
*
*  3.2 Case of Leap-frog - correct velocities due to barostat
*
      HTSTEP=0.5*TSTEP
      if(IAL.eq.3.and.LNVT)then
*   3.2.1    Separate thermostat on each molecule type
*            (velocity scaled before thermostat parameter correction)
        if(LSCTD)then
          do I=1,NSTOT
	    ITYP = ITYPE(I)
            SCV= 1.-SCM(ITYP)*TSTEP
	    if(LMOVE(ITYP))then
              VX(I) = VX(I)*SCV
              VY(I) = VY(I)*SCV
              VZ(I) = VZ(I)*SCV
            end if
          end do
*   3.2.2  Common thermostat
        else
          SCV = 1.-SC*TSTEP
          do I=1,NSTOT
	    ITYP = ITYPE(I)
	    if(LMOVE(ITYP))then 
              VX(I) = VX(I)*SCV
              VY(I) = VY(I)*SCV
              VZ(I) = VZ(I)*SCV
            end if
          end do
        end if
      end if
*
*  3.3 Recalculate temperatures and pressures
*  ---------------
*  3.3.1 Calculate temperature    
      call GETEMP 
*  3.3.2 Calculate pressure  
C        True only for MAST node
C        This collect virial for "atomic" pressure
      VIRSUM	= VIR1+VIR2+VIRB+PELS1+PELS2+SPE+VIRD+
     +            VIRA(1)+VIRA(2)+VIRA(3)
C        This collect virial for "molecular" pressure
      WIRSUM	= VIR1+VIR2+PELS1+PELS2+SPE+VIRD+WIRS+WIRSS
CD      write(*,'(5e13.5)')VIR1,VIR2,PELS1,PELS2,SPE,VIRD,WIRS,WIRSS
C        Kinetic (ideal gas) contributions to pressure
      TRID	= 2.*TKE*UNITP/(3.*VOL)
      TRIDM	= NOP*BOLTZ*TEMP*1.d-5/(VOL*UNITL**3)
C        Pressure in atm and its projections
      TRYCK	= VIRSUM*UNITP/(3.*VOL)+TRID
      TRYCKM	= WIRSUM*UNITP/(3.*VOL)+TRIDM 
C        Projections are calculated for "atomic" pressure
      TRYCKX	= (VIRX+VIRFX+VIRD/3.+VIRA(1))*UNITP/VOL+TRID
      TRYCKY	= (VIRY+VIRFY+VIRD/3.+VIRA(2))*UNITP/VOL+TRID
      TRYCKZ	= (VIRZ+VIRFZ+VIRD/3.+VIRA(3))*UNITP/VOL+TRID
CD	if(TASKID.eq.MAST)write(*,'(5f13.2)')
CD     +TRYCK,TRYCKX,TRYCKY,TRYCKZ,(TRYCKX+TRYCKY+TRYCKZ)/3.
C
C     Now temperature TEMP and pressure TRYCK are defined for the 
C     current configuration, as well as their components
*
*   3.3.3 Calculate deviation from thermostat T and P
      DKE         = TEMP/TRTEMP-1.d0               !  total temperature 
      do I=1,NTYPES
CD         write(*,*)' temp ',I,TEMPR(I),TASKID
         DTE(I)   = TEMPR(I)/TRTEMP-1.d0           !  for each species
      end do
C    in the case of constrained dynamics, pressure is defined by 
C    "molecule" algorithm - it corresponds scaling of molecular COM
C    for flexible molecules, scaling of atom positions is employed,
C    which corresponds to the "atomic" pressure 
C    Exception: case of separate pressure control in each direction 
C    (LSEP=.true.). Then pressure in each direction is determined from
C    "atomic" pressure. 
      if(LSHEJK)then
	 DPE(1) = TRYCKM-TRPRES
      else
         DPE(1) = TRYCK -TRPRES
      end if
      DPE(2) = TRYCKX-TRPRES
      DPE(3) = TRYCKY-TRPRES
      if(LHEX)then
        DPE(2)=0.5*(DPE(2)+DPE(3))