mdstepe.f

分子动力学模拟程序

        DPE(3)=DPE(2)
      end if
      DPE(4) = TRYCKZ-TRPRES
C    This broadcast pressure to all nodes
*
*   3.4  First integration of Nose-Hoover equations  (t->t+1/2)
*   -----------------------------------------------------------
      if(IAL.eq.1)then
        if(LNVT)then
*   3.4.1 Correction ksi due to thermostat
C         coefficients DQT,DQP and RTP were defined in main.f, p.6.2
C
C    Nose eqn is:    dP/dt   -> F - P*ksi/Q
C                    dksi/dt -> 2*Ekin - Nf*kT
C    Here ksi/Q = SC
C             Q = Nf*kT*tau**2  ; tau -> QT (i.u.)
C                                 Q/(Nf*kT) -> 1./DQT 
C    So     d(SC)/DT = (2*Ekin/Nf*kT-1)*DQT = DKE*DQT
C
C    This is separate thermostat for each species 
          if(LSCTD)then
            SCA = 0.
            do I=1,NTYPES
              SCM(I) = SCM(I) + DTE(I)*DQT*HTSTEP     !  DQT -thermostat mass
              SCA = SCA + SCM(I)*EKIN(I) 
            end do
            SC = SCA/TKE
C   common thermostat
          else
 	    SC = SC + DKE*DQT*HTSTEP
          end if
          if(LNPT)then
*   3.4.2 Barostat corrections (NPT)
*   -------------------------------
*   3.4,2,1 Correction ksi due to barostat
C
C   Additional correction to ksi:
C       d(ksi)=ita**2/W-kT
C       ita/W -> SCL
C     so d(ksi)  ->   SCL**2*(W/Q) - kT/Q 
C                       
            if(LSCTD)then      
              SCA = 0.
              do I=1,NTYPES
                SCM(I) = SCM(I) + (SCL**2*RTP-DQT/FNST)*HTSTEP
                SCA = SCA + SCM(I)*EKIN(I) 
              end do
              SC = SCA/TKE
            else
              SC = SC + (SCL**2*RTP-DQT/FNST)*HTSTEP    !  RTP=DQT/DQP
            end if
*   3.4.2.2 Correction ita 
C         TKE is kinetic energy
            if(LSEP)then
              SCLX = SCLX+DQP*(VOL*DPE(2)+2.*TKE/FNST-SCLX*SC)*HTSTEP
              SCLY = SCLY+DQP*(VOL*DPE(3)+2.*TKE/FNST-SCLY*SC)*HTSTEP
              if(LHEX)then
                 SCLX=0.5*(SCLX+SCLY)
                 SCLY=SCLX
              end if
              SCLZ = SCLZ+DQP*(VOL*DPE(4)+2.*TKE/FNST-SCLZ*SC)*HTSTEP
              SCL = SCLX+SCLY+SCLZ   ! trace(P)
            else
              SCL = SCL + DQP*(3.*VOL*DPE(1)+6.*TKE/FNST-SCL*SC)*HTSTEP
              SCLX=SCL
              SCLY=SCL
              SCLZ=SCL
C    Control fluctuations
	      if(abs(SCL*HTSTEP).ge.0.1)then
	        write(*,*)' too strong fluctuations in NPT algorithm'
	        write(*,*)' scaling factor ',SCL*HTSTEP
   	        if(IPRINT.ge.7)then
	          write(*,*)'------------------->'
	          write(*,*)SCL,TRYCK,TRYCKM
	          write(*,*)VIR1*UNITP/(3.*VOL),VIR2*UNITP/(3.*VOL),
     + PELS1*UNITP/(3.*VOL),PELS2*UNITP/(3.*VOL),VIRB*UNITP/(3.*VOL)
	          write(*,*)'<-------------------'
	        end if
                ISKK=ISKK+1
                if(ISKK.gt.100)then 
	write(*,*)' !!! repeated failure of NPT-algorithm' 
	write(*,*)' !!! restart program with constant volume'
	write(*,*)' !!! or increase thermostat parameter for pressure'
                  stop
                end if
              end if    ! if(LSEP   
            end if
*   3.4.2.3 Calculate scaling coefficients  
            SCVX = 1.-((1.-3./FNST)*SCL+SC)*HTSTEP
            SCVY = SCVX
            SCVZ = SCVX
            DRCX = TSTEP*SCLX
            DRCY = TSTEP*SCLY
            DRCZ = TSTEP*SCLZ
            DRC  = (DRCX+DRCY+DRCZ)/3.
            DRCV = DRC*HTSTEP           !  These are second order corrections 
            DRCX = 1.+DRCX + 0.5*DRCX**2 !  and (may be?) ommitted
            DRCY = 1.+DRCY + 0.5*DRCY**2 !  and (may be?) ommitted
            DRCZ = 1.+DRCZ + 0.5*DRCZ**2 !  and (may be?) ommitted
*
*   3.4.2.4 Correct coordinates and velocities - flexible molecules
            do I=1,NSTOT
	        ITYP = ITYPE(I)
	        if(LMOVE(ITYP))then 
                  VX(I) = VX(I)*SCVX
                  VY(I) = VY(I)*SCVY
                  VZ(I) = VZ(I)*SCVZ
                  SX(I) = SX(I)*DRCX + VX(I)*DRCV*MASSDI(I)
                  SY(I) = SY(I)*DRCY + VY(I)*DRCV*MASSDI(I)
                  SZ(I) = SZ(I)*DRCZ + VZ(I)*DRCV*MASSDI(I)
                end if
            end do
*   3.4.2.5  Scale simulation box
            call RECLEN(DRCX,DRCY,DRCZ)
         else      !  .not.LNPT
*  3.4.3 Velocities corrections in the NVT ensemble 
           call GETEMP
           if(LSCTD)then
             do I=1,NSTOT
	       ITYP = ITYPE(I)
               SCV = 1.-SCM(ITYP)*HTSTEP
	       if(LMOVE(ITYP))then
                 VX(I) = VX(I)*SCV
                 VY(I) = VY(I)*SCV
                 VZ(I) = VZ(I)*SCV
               end if
             end do
           else
             SCV = 1.-SC*HTSTEP
             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   ! if(LSCTD)
         end if     ! if(LNPT)
*   3.4.4  Absorbed energy
         EABS = EABS - (SC+SCL*(1.-3./FNST))*TKE*TSTEP
       end if       ! if(LNVT)
*
*   3.5  Second integration of Nose-Hoover equations (t+1/2->t+1)
*   -------------------------------------------------------------
      else if(IAL.eq.2)then
        if(LNVT)then
*   3.5.1 Correct velocities due to thermostat
          if(LSCTD)then
            do I=1,NSTOT
	      ITYP = ITYPE(I)
              SCV = 1.-SCM(ITYP)*HTSTEP
	      if(LMOVE(ITYP))then
                VX(I) = VX(I)*SCV
                VY(I) = VY(I)*SCV
                VZ(I) = VZ(I)*SCV
              end if
            end do
          else
            SCV = 1.-SC*HTSTEP
            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
*   3.5.2 Correct velocities due to barostat
          if(LNPT)then
            SCVX = 1.-(1.-3./FNOP)*SCL*HTSTEP
            SCVY = SCVX
            SCVZ = SCVX
            do I=1,NSTOT
	      ITYP = ITYPE(I)
	      if(LMOVE(ITYP))then 
                VX(I) = VX(I)*SCVX
                VY(I) = VY(I)*SCVY
                VZ(I) = VZ(I)*SCVZ
              end if
            end do
*   3.5.3  Recalculate temperature
            call GETEMP 
*   3.5.4  Correct ita
            if(LSEP)then
              SCLX = SCLX+DQP*(VOL*DPE(2)+2.*TKE/FNST-SCLX*SC)*HTSTEP
              SCLY = SCLY+DQP*(VOL*DPE(3)+2.*TKE/FNST-SCLY*SC)*HTSTEP
              if(LHEX)then
                 SCLX=0.5*(SCLX+SCLY)
                 SCLY=SCLX
              end if
              SCLZ = SCLZ+DQP*(VOL*DPE(4)+2.*TKE/FNST-SCLZ*SC)*HTSTEP
              SCL = SCLX+SCLY+SCLZ   ! trace(P)
            else
              SCL = SCL + DQP*(3.*VOL*DPE(1)+6.*TKE/FNST-SCL*SC)*HTSTEP
              SCLX=SCL
              SCLY=SCL
              SCLZ=SCL
C    Control fluctuations
	      if(abs(SCL*HTSTEP).ge.0.1)then
	        write(*,*)' too strong fluctuations in NPT algorithm'
	        write(*,*)' scaling factor ',SCL*HTSTEP
   	        if(IPRINT.ge.7)then
	          write(*,*)'------------------->'
	          write(*,*)SCL,TRYCK,TRYCKM
	          write(*,*)VIR1*UNITP/(3.*VOL),VIR2*UNITP/(3.*VOL),
     + PELS1*UNITP/(3.*VOL),PELS2*UNITP/(3.*VOL),VIRB*UNITP/(3.*VOL)
	          write(*,*)'<-------------------'
	        end if
                ISKK=ISKK+1
                if(ISKK.gt.100)then 
	write(*,*)' !!! repeated failure of NPT-algorithm' 
	write(*,*)' !!! restart program with constant volume'
	write(*,*)' !!! or increase thermostat parameter for pressure'
                  stop
                end if
              end if    ! if(LSEP   
            end if
*   3.5.5  Correct ksi due to barostat
            SC = SC + (SCL**2*RTP-DQT/FNST)*HTSTEP
            if(LSCTD)then
              do I=1,NTYPES
                SCM(I) = SCM(I) + (SCL**2*RTP-DQT/FNST)*HTSTEP
              end do
            end if
          else      !  .not.LNPT
*   recalculate temperature
            call GETEMP
          end if    ! if(LNPT
*   3.5.6 Correction ksi due to thermostat
          if(LSCTD)then
            SCA = 0.
            do I=1,NTYPES
              DTE(I)   = TEMPR(I)/TRTEMP-1.d0          
              SCM(I) = SCM(I) + DTE(I)*DQT*HTSTEP
              SCA = SCA + SCM(I)*EKIN(I) 
            end do
            SC = SCA/TKE
          else 
            DKE  = TEMP/TRTEMP-1.d0
	    SC = SC + DKE*DQT*HTSTEP
          end if
*   3.5.7  Absorbed energy
         EABS = EABS - (SC+SCL*(1.-3./FNST))*TKE*TSTEP
        end if ! if (LNVT
*
*   3.6  Case of single-step SHAKE algorithm
*   ----------------------------------------
      else if(IAL.eq.3)then
*   3.6.1  Correct ksi(t+1) from velocities and temperature at t+1/2
        if(LNVT)then
	  SC = SC + DKE*DQT*TSTEP      !  Ksi(t+1)
C    This is separate thermostat for each species (may be not work)
          if(LSCTD)then
            SCA = 0.
            do I=1,NTYPES
              SCM(I) = SCM(I) + DTE(I)*DQT*TSTEP     !  DQT -thermostat mass
              SCA = SCA + SCM(I)*EKIN(I) 
            end do
            SC = SCA/TKE
          end if
          if(LNPT)then
*   3.6.2 Barostat corrections (NPT)
*   -------------------------------
*   3.6.2.1 Correction ita  t-1/2  -> t+1/2 
C         TKE is kinetic energy
            if(LSEP)then
              SCLX = SCLX+DQP*(VOL*DPE(2)+2.*TKE/FNST-SCLX*SC)*TSTEP
              SCLY = SCLY+DQP*(VOL*DPE(3)+2.*TKE/FNST-SCLY*SC)*TSTEP
              if(LHEX)then
                 SCLX=0.5*(SCLX+SCLY)
                 SCLY=SCLX
              end if
              SCLZ = SCLZ+DQP*(VOL*DPE(4)+2.*TKE/FNST-SCLZ*SC)*TSTEP
              SCL = SCLX+SCLY+SCLZ   ! trace(P)
            else
              SCL = SCL + DQP*(3.*VOL*DPE(1)+6.*TKE/FNST-SCL*SC)*TSTEP
              SCLX=SCL
              SCLY=SCL
              SCLZ=SCL
*   3.6.2.2    Control fluctuations
	      if(abs(SCL*TSTEP).ge.0.1)then
	        write(*,*)' too strong fluctuations in NPT algorithm'
	        write(*,*)' scaling factor ',SCL*TSTEP
   	        if(IPRINT.ge.7)then
	          write(*,*)'------------------->'
	          write(*,*)SCL,TRYCK,TRYCKM
	          write(*,*)VIR1*UNITP/(3.*VOL),VIR2*UNITP/(3.*VOL),
     + PELS1*UNITP/(3.*VOL),PELS2*UNITP/(3.*VOL),VIRB*UNITP/(3.*VOL)
	          write(*,*)'<-------------------'
	        end if
                ISKK=ISKK+1
                if(ISKK.gt.100)then 
	write(*,*)' !!! repeated failure of NPT-algorithm' 
	write(*,*)' !!! restart program with constant volume'
	write(*,*)' !!! or increase thermostat parameter for pressure'
                  stop
                end if
              end if    ! if(LSEP   
            end if
*  this is for correction of velocities
            SCV = ((1.-3./FNST)*SCL)*TSTEP
*   3.6.2.3 Calculate scaling coefficients  
            DRCX = TSTEP*SCLX       ! eta*dt
            DRCY = TSTEP*SCLY
            DRCZ = TSTEP*SCLZ
*   3.6.2.4 Correct coordinates and velosities 
*   3.6.2.4.1   Calculate centre of mass and molecular momenta
C    Cycle over molecules
	    do IMOL = 1,NOP
              ITYP = ITM(IMOL)
              SUMM        = SUMMAS (ITYP)
	      if(LMOVE(ITYP))then 
		NSBEG = ISADR(ITYP)+1
		NSEND = ISADR(ITYP+1)
		XC       = 0.D0
      	        YC       = 0.D0
                ZC       = 0.D0
	        PXC	= 0.
                PYC	= 0.
                PZC	= 0.
C   Cycle over atoms within the molecule                  
		IBEG=ISADDR(ITYP)+1+NSITS(ITYP)*(IMOL-IADDR(ITYP)-1)
		IEND=ISADDR(ITYP)+NSITS(ITYP)*(IMOL-IADDR(ITYP))
                DO I       = IBEG,IEND
                  IS	     = NSITE(I)
                  XC      = XC+MASS(IS)*SX(I)
                  YC      = YC+MASS(IS)*SY(I)
                  ZC      = ZC+MASS(IS)*SZ(I)
C   Molecular momenta
                  PXC     = PXC+VX(I)
                  PYC     = PYC+VY(I)
                  PZC     = PZC+VZ(I)
                END DO! OF N
                XC       = XC/SUMM    
      	        YC       = YC/SUMM
      	        ZC       = ZC/SUMM
*   3.6.2.4.2  Correction to COM momenta
                DVX	= PXC*SCV
                DVY	= PYC*SCV
                DVZ	= PZC*SCV 
*   3.6.2.4.3  COM displacements
                DXC = XC*DRCX
                DYC = YC*DRCY
                DZC = ZC*DRCZ
*   3.6.2.4.4  Correct atom velocities
C    DVX - molecular momentum correction
C    DVX/SUMM - COM velocity correction
C    MASS(IS)*DVX/SUMM - atom momentum correction
      	        DO I       = IBEG,IEND
      	          IS	= NSITE(I)
                  FAC	= MASS(IS)/SUMM
                  VX(I)	= VX(I)-DVX*FAC
                  VY(I)	= VY(I)-DVY*FAC
                  VZ(I)	= VZ(I)-DVZ*FAC
*   3.6.2.4.5  Correct atom positions
                  SX(I) = SX(I) + DXC 
                  SY(I) = SY(I) + DYC 
                  SZ(I) = SZ(I) + DZC 
                end do
	      end if
            end do ! of ITYP
*   3.6.2.5  Scale simulation box
            DRCX=DRCX+1.
            DRCY=DRCY+1.
            DRCZ=DRCZ+1.
            call RECLEN(DRCX,DRCY,DRCZ)
*   3,6,2,6 Correction ksi due to barostat
C
C   Additional correction to ksi:
C       d(ksi)=ita**2/W-kT
C       ita/W -> SCL
C     so d(ksi)  ->   SCL**2*(W/Q) - kT/Q 
C                       
            SC = SC + (SCL**2*RTP-DQT/FNST)*TSTEP    !  RTP=DQT/DQP
            if(LSCTD)then    
              SCA = 0.
              do I=1,NTYPES
                SCM(I) = SCM(I) + (SCL**2*RTP-DQT/FNST)*TSTEP
                SCA = SCA + SCM(I)*EKIN(I) 
              end do
              SC = SCA/TKE
            end if
          end if     ! if(LNPT)
*   3.6.4  Absorbed energy
          EABS = EABS - (SC+SCL*(1.-3./FNST))*TKE*TSTEP
        end if       ! if(LNVT)
      end if   ! if(IAL
*
*  3.7 Forcible adjusting of velocities to given temperature
*  ---------------------------------------------------------
      IF(LSCLT.and.(IAL.eq.2.or.IAL.eq.3)) THEN 
        if(LSCTD)then
          DO ITYP     = 1,NTYPES
            IF(DABS(TRTEMP-TEMPR(ITYP)).GT.TDELT.and.LMOVE(ITYP))THEN
              SCV = DSQRT(TRTEMP/TEMPR(ITYP))
	      IBEG=ISADDR(ITYP)+1
	      IEND=ISADDR(ITYP+1) 
	      do I=IBEG,IEND
                VX(I)       = SCV*VX(I)
                VY(I)       = SCV*VY(I)
                VZ(I)       = SCV*VZ(I) 
	      end do  
              if(IPRINT.ge.5)write(*,*)
     +      'velocities scaled by ',SCV,' for type',ITYP
              NRTSC(ITYP) =  NRTSC(ITYP)+1
	    end if
	  end do
        else
          IF(DABS(TRTEMP-TEMP).GT.TDELT)THEN
            SCV = dsqrt(TRTEMP/TEMP)
            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
            if(IPRINT.ge.5)write(*,*)
     +      'velocities scaled by ',SCV
            NRTSC(ITYP) =  NRTSC(ITYP)+1
          end if
        end if
        call GETEMP
      END IF!(LSC...
      TKE	= TKE/FNOP
*
      RETURN
      END