mdee.f

分子动力学模拟程序

        DQP	= UNITE*UNITT**2/(UNITP*FNOP*BOLTZ*TRTEMP*(QPR*1.d-15)**2)
        write(*,'(a,f10.3,a)')
     +  '*** Constant-pressure MD: P = ',TPRES,' atm'
        PRINT *,'   NOSE-NPT RELAXATION TIME     ',QPR,' fs'
        if(LSEP)write(*,*)
     +  ' Separate pressure controll in different directions'
      else
        LSEP=.false.
      END IF!(LNVT...
        SC	 = 0.D0
        SCL	= 0.D0
        SCLX = 0.
        SCLY = 0.
        SCLZ = 0.        
        do I=1,NTYPES
           SCM(I)=0.
        end do
*
Continue an old run:
*
      IF(LRST) THEN
*
         BOXLN=BOXL
         BOYLN=BOYL
         BOZLN=BOZL         
        CALL RESTRT(1)
        call RECLEN(ODIN,ODIN,ODIN)
      call RECINT
*     
	if(LECH)then
	  do I=1,NE
	    EE(I)=EE(I)*NDEL
	  end do
	end if 
      CALL GETCOM
      write(*,'(a,f12.4,a2)')
     +' The system was simulated ',TIM*1.d12,'ps'
      write(*,'(a,f12.4,a2)')
     +' Averages was collected   ',TIMA*1.d12,'ps'
*
      ELSE
*
*Prepare a new run:
*
*Generate start coordinates for C.O.M. points:
*
      NHIST	= 1
      TIM	= 0.
      timest	= 0.
	ICHE=0
	ICHU=0
	NTRKF=0           ! num. of trace file
	NTRKZ=0           ! num. of points in the trace file
      IF(NOP.EQ.1) THEN
      X(1)   = 0.D0
      Y(1)   = 0.D0
      Z(1)   = 0.D0
*
      ELSE
*
	  if(ICHECK.eq.0)then
*  7.2.2.5  read atom coord. from input file
	    open(unit=13,file=FINP,status='old',err=298)
	    go to 299
C   something goes wrong:
C   can not read input file...         
 297	      if(TASKID.eq.MAST)then
              write(*,*)'too few atoms in the input file'
	        IE=99
	        AUX=TAKESTR(13,IE)
	        write(*,*)' fatal processor error'
	      end if
	      stop       
C   can not open input file                    
 298	      if(TASKID.eq.MAST)
     +write(*,*)' input file ',FINP,' not found (ICHECK=0)'
	      stop
C    reading input in XMOL format
 299          read(13,*)NNN
              if(NNN.ne.NSTOT)write(*,*)
     + '!!! Obs: number of atoms in coordinate file differ from input:',
     +  NNN,' and ',NSTOT
              read(13,*)
              do I=1,NSTOT                                
		AUX=TAKESTR(13,IE)
	        read(AUX,*,err=297,end=297)CDUM,SX(I),SY(I),SZ(I) 
	      end do          
	      if(TASKID.eq.MAST)write(*,*)
     +'read initial coordinetes from file ',FINP
	      close(13)  
              call GETCOM
*  7.2.2.6  read molecular center-of-mass coord. from input file
	    else if(ICHECK.eq.-1)then
	      open(unit=13,file=FINP,status='old',err=398)
	      go to 399         
C   something goes wrong:
C   can not read input file...         
 397	      if(TASKID.eq.MAST)then
              write(*,*)'too few molecules in the input file'
	        IE=99
	        AUX=TAKESTR(13,IE)
	        write(*,*)' 1.e1000 thanks to your system administrator'
	      end if
	      stop
C   can not open input file                    
 398	      if(TASKID.eq.MAST)
     +write(*,*)' input file ',FINP,' not found (ICHECK=-1)'
	      stop
C    reading...
 399          do I=1,NOP                                
		AUX=TAKESTR(13,IE)
	        read(AUX,*,err=397,end=397)X(I),Y(I),Z(I)
		call PBC(X(I),Y(I),Z(I))
	      end do          
	      if(TASKID.eq.MAST)write(*,*)
     +'read initial COM from file ',FINP
	      close(13)  
            else IF(ICHECK.eq.1) THEN
              CALL FCC(X,Y,Z,BOXL,NOP)
            ELSE
              CALL CUBIC(X,Y,Z,FX,FY,FZ,BOXL,NOP)
            END IF!(ICHE...
*
      END IF!(NOP...
*
*Initialize velocities
*
      TMPRED     = EPSIL(1)*504.D0/TRTEMP
      CALL COMVEL
*
      PRINT *,'*** RANDOM VELOCITIES FOR   ',NOP,'  PARTICLES'
*
      VSQ         = 0.0
      DO N        = 1,NSTOT
      VSQ         = VSQ+(VX(N)**2+VY(N)**2+VZ(N)**2)*MASSDI(N)
      END DO! OF N
*
      TEMP        = 1.50*VSQ*ENERF*TEMPF/FNST
      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! OF N
*
      VSQ         = 0.0
      DO N        = 1,NSTOT
      VSQ         = VSQ+(VX(N)**2+VY(N)**2+VZ(N)**2)*MASSDI(N)
      END DO! OF N
*
      TEMP        = 1.5*VSQ*ENERF*TEMPF/FNST
*
      PRINT *,'*** INITIAL TEMPERATURE CALCULATED',TEMP
      PRINT *,'*** VELOCITIES INITIATED'
      call GETEMP
      PRINT *
*
*Orientations of the molecules:
*
      if(ICHECK.ne.0)then
      N              = 0
      I              = 0
      DO ITYP        = 1,NTYPES
      IBES           = ISADR (ITYP)+1
      IENS           = ISADR (ITYP +1)
      DO J           = 1,NSPEC(ITYP)
      I              = I+1
      CALL RANQ(QP)
      DO IS          = IBES,IENS
      N              = N+1
      RXI            = R(IS,1)
      RYI            = R(IS,2)
      RZI            = R(IS,3)
      CALL ROTATE  (QP,RXI,RYI,RZI)
      SX(N)          = RXI+X(I)
      SY(N)          = RYI+Y(I)
      SZ(N)          = RZI+Z(I)
      WX(N)          = RXI
      WY(N)          = RYI
      WZ(N)          = RZI
      GX(N)          = 0.D0
      GY(N)          = 0.D0
      GZ(N)          = 0.D0
      HX(N)          = 0.D0
      HY(N)          = 0.D0
      HZ(N)          = 0.D0
      OX(N)          = SX(N)
      OY(N)          = SY(N)
      OZ(N)          = SZ(N)
      END DO! OD IS
      END DO! OF I
      END DO! OF ITYP
*
*
      PRINT "('*** MOLECULAR ORIENTATIONS INITIALIZED')" 
      end if
      NSTTOT	= 0      
      call RECLEN(ODIN,ODIN,ODIN)	    
	call RECINT
	  do I=1,NE
	    EE(I)=EE(I)*NDEL
	  end do
*
      END IF!(LRST...
* Parameter for Ewald sum
	if(RCUTT.le.0.5*BOXL)then
        ALPHAD	= ALPHA/RCUTT
	else
	  ALPHAD	= 2.d0*ALPHA/BOXL 
	end if 
*
*Zero average accumulators
*     
	if(L0AVS)then
	  NAVT	= 0
	  write(*,*)' New averaging'
          NSTEPA	= 0
          TIMA	= 0.
	  call ZEROAV 
          NSTEPA=0
          DO I        = 1,NRQS
	    do J=1,NE
              AV(I,J)       = 0.D0
              AW(I,J)       = 0.D0
	    end do
          END DO! OF I
          do I=1,NE
            do J=1,NE
              IWALK(I,J)=0
            end do
          end do
	end if
	TRPRES=TPRES
	PRINT "(/'*** DENSITY           ',D12.6,' g/cm**3 ' )",RHO
      write(*,'(a,3f12.4)')'*** Box sizes:    ',BOXL,BOYL,BOZL
      PRINT "( '*** BOX VOLUME        ',D12.6,' A**3    ' )",VOL
      PRINT "( '*** BOX LENGTH        ',D12.6,' A       ' )",BOXL
      PRINT "( '*** CUT-OFF           ',D12.6,' A       ' )",RCUTT
	PRINT "( '*** Closest R         ',D12.6,' A       ' )",SHORT
	PRINT "( '*** ALPHA (Ew.sum)    ',D12.6)",ALPHA  
	PRINT "( '*** Fexp  (Ew.sum)    ',D12.6)",FEXP
      PRINT *       
*
      IF(LGR)   CALL RDFINP
*  prepare list of bound atoms
      CALL BNDLST
*  change/check  temperature/density
      call CHTERM(BOXLN,BOYLN,BOZLN)
      call RECLEN(ODIN,ODIN,ODIN)
      call TESTCONF
      IF(LCHECK)then
	  if(LRST)then
	    call GETAVR(3)
	    IF(LGR.and.LGRST)   CALL RDFOUT(1)
        else  
	    WRITE(*,*)' STOP after cheking the input' 
	  end if
	  if(LXMOL)call XMOLDUMP
	  go to 9999
      end if
	PRINT "('*** MD STEPS ABOUT TO START FOR',I8,' STEPS:')",NSTEPS
      PRINT * 
*
      if(IPRINT.ge.5)WRITE(6,900)
* A MD step 
      IF(LSHEJK) THEN
        CALL SULTAN
      ELSE
        CALL DOUBLE
      END IF
*
      if(LXMOL)call XMOLDUMP
      if(IHIST.gt.NHIST)IHIST=NHIST
      call GETAVR(3)
      IF(LGR.and.LGRST)   CALL RDFOUT(1)
*
  900 FORMAT(/1X,'  STEP     POTE      PE(int)   E(el)'
     X,'   Etot       TEMP     PRES      Dens '/)
*
      IF(.NOT.LNVT) THEN
      PRINT *,'... TEMP SCALED:'
      DO ITYP = 1,NTYPES
      PRINT *,'...',NRTSC(ITYP),'  TIMES IN'
     X,NSTEPS,'  STEPS - FOR TYPE: ',ITYP
      END DO! OF ITYP
      END IF
 9999 stop
*
      END
*