mdee.f

分子动力学模拟程序

      PROGRAM MDEE
*
      include "mdee.h"
*-----------------------------------------------
      DIMENSION  QP(4)
      character*4 CDUM
*
      NAMELIST /INPUT/ NTYPES,NSPEC,NSITS,FNAME,TPRES
     X         ,TRTEMP,RHO,DT,RCUT,TDELT,FTSCL,QT,QPR
     X         ,BOXL,BOYL,BOZL,BDTOL,TOLER,C14LJ,C14EL
     X         ,NSTEPS,NFREQ,IPOT,NAVER,IPRINT,NAMOL
     X         ,NDUMP,IAVER,ICHECK,IHIST,ALPHA,FEXP
     X         ,LRST,LDMP,LSCLT,LSCFT,LNVT,LNPT,LCHECK,L14NB
     X         ,L0AVS,L0VS,LSHEJK,LCHP,LCHT,LXMOL,LECH,L15NB
     X         ,SHORT,ICHNB,LGR,LGDMP,LGRST,ALLSTS,RDFCUT,PATHDB
     X         ,LCFDMP,LCFRST,NSTEG,JUMP,N1,N2,MAX
     +         ,ME0,NE,IEE,NCEE,PLOW,LSEP,LSCTD
*  
      character*20 FNM, FINP*24,PATHDB*64,DATE*8
      character*128 TAKESTR,AUX
      data ODIN/1.d0/
      label = "MDEE_v4.4   "
	DATE  = '30.11.05'
      PRINT *,'*** MOLECULAR DYNAMICS',
     X' PROGRAM FOR MIXTURES OF RIGID or FLEXIBLE MOLECULES'
      PRINT *,'*** ',
     X '         with removing/inserting of molecules '
	write(*,*)'     (expanded ensemble method)'
      PRINT *
	write(*,*)'Version ',label,'           from ',DATE
	write(*,*)
	time0	= cputime(0.d0)
*
*     Default parameters
*
*Physical parameters (default) for a many-particle system:
*
      RHO        = 1.00              ! density   (if BOXL=0)
      TRTEMP     = 300.0             ! temperature
      TPRES      = 1.                ! pressure (atm)
      NOP        = NPART             ! number of particles
      NFREQ      = 10                ! num. of short time steps in DT
      NAVER      = 100               ! num. of steps for averaging 
      DT         = 1.D-15            ! time step
      RCUT       = 10.D0             ! cut-off
      TDELT      = 50.D0             ! ??
      RDFCUT     = 10.D0             ! cut off for RDF  
      ALPHA	= 2.5                  ! Ewald: erfc(alpha)=0
      FEXP	= 8.                   !        exp(-fexp)=0
      NUMTASK=1
      TASKID=0
      LHEX=.false.
      LOCT=.false.
      FTSCL      = 0.001                ! cutoff large forces 
      QT	 = 10.D0
      QPR	 = 300.d0
      JUMP       = 1
      BOXL       = 0.D0
      BOYL       = 0.D0
      BOZL       = 0.D0
      SHORT	= 5.                   ! R for closest neigbours   
      TOLER      = 1.D-4
      BDTOL      = 1.D-4 
	MAX	= MAXS
      ICHECK	= 1
      LMNBL=.false.
	IHIST	= 1
	ICHNB	= 10    ! cheking neigbours   list
	NCEE	= 1     ! num. of steps followed by change of subensemble
*
*Some control parameters
*
      FNAME	 = "md"
      LCHECK     = .FALSE.              ! true - ony input checking
      LNVT       = .FALSE.
      LNPT       = .FALSE.
      LSCLT      = .FALSE.
      LSCFT      = .FALSE.
      L0AVS      = .FALSE.
      L0VS       = .FALSE.
      LRST       = .FALSE.
      LDMP       = .FALSE.
      LCHL       = .FALSE.
      LCHP       = .FALSE.
      LCHT       = .FALSE.
      LECH	= .false.
      LSEP      = .false.
      LSCTD=.false.
      MSTEP      =  0
*     
      DO I       = 1,NTPS
        ISHEJK(I)  = 1
        NRB(I)     = 0
        C14LJ(I)   = 1.D0
        C14EL(I)   = 1.D0
        L14NB(I)=.true.
        L15NB(I)=.true.
        LMOVE(I)=.true.
      END DO
*
*Read in new values
*
      READ(*,INPUT)
      if(.not.LRST)then
	  L0AVS=.true.
	  ME=ME0
      end if
      do IE=1,NE
	read(*,*)EC(IE),EE(IE) 
	if(IPRINT.ge.7)write(*,*)EC(IE),EE(IE)
*  some historical reasons...
        EE(IE)=-EE(IE)
      end do    
      if(LCHECK)L0AVS=.false.
      if(MAX.gt.MAXS)stop "increase MAXS"
      FTSCL=1.d0/FTSCL
*
*   file names
*
      IC=1
      do I=1,20
        if(FNAME(I:I).ne.' ')IC=I
      end do
      IC1=IC+1
      IC2=IC+4
	ICFN=IC+2
      FNM='                    '
	write(FDUMP,'(24(1h  ))')
	write(FLST,'(24(1h  ))')
	write(FRDF,'(24(1h  ))')
	write(FTRK,'(24(1h  ))')
      write(FBIO,'(24(1h  ))')
      write(FINP,'(24(1h  ))')
      FNM(1:IC)=FNAME(1:IC)
      FDUMP(1:IC)=FNM
      FDUMP(IC1:IC2)='.dmp'
      FLST(1:IC)=FNM
      FLST(IC1:IC2)='.lst'
      FRDF(1:IC)=FNM
      FRDF(IC1:IC2)='.rdf'
      FBIO(1:IC)=FNM
      FBIO(IC1:IC2+1)='.xmol'
      FINP(1:IC)=FNM
      FINP(IC1:IC2)='.inp'
	CLN=BOXL
*
*---------------- DEFINE MOLECULE TYPES -------------------------------
*
      SUMM         = 0.D0
      NRBB         = 0
      NRAA         = 0
      NRTT         = 0
      NRII         = 0
      NX           = 0
      IADB(1)      = 1
      IADA(1)      = 1
      IADT(1)      = 1
      TOTMAS       = 0.D0
	NEFAC=0
*
*  Input of molecules from data base
        NTYP=0
*     
	do ITYP=1,NTYPES
        CALL READMOL (SUMM,NRBB,NRAA,NRTT,NRII,NX,NTYP,PATHDB,
     +  NAMOL(ITYP))
        TOTMAS       = TOTMAS+SUMM*NSPEC(ITYP)
        SUMMAS(NTYP) = SUMM
	  NAME(ITYP)=NAMOL(ITYP)(1:6)
        if(IEE(ITYP).eq.0)NEFAC=NEFAC+NSPEC(ITYP)
      end do
*     
      IF(NX  .GT.NS  ) STOP '!!! INCREASE NS   !'
      IF(NTYPES.GT.NTPS) STOP '!!! INCREASE NTPS !'
      IF(NRBB.GT.NB  ) STOP '!!! INCREASE NB   !'
      IF(NRAA.GT.NA  ) STOP '!!! INCREASE NA   !'
      IF(NRTT.GT.NT  ) STOP '!!! INCREASE NT   !'
      IF(NRII.GT.NI  ) STOP '!!! INCREASE NI   !'
      IF(TOTMAS.EQ.0 ) STOP '!!! NO PARTICLES WITH MASSES DEFINED !'
*
      IF(LSHEJK) THEN
      NFREQ = 1
      PRINT "(/'*** CONSTRAINED DYNAMICS WILL BE USED! ***'/)"
      END IF
	NRBT=0.
	do ITYP=1,NTYPES
	  NRBT=NRBT+NRB(ITYP)*NSPEC(ITYP)
	end do
*
      NOP     = 0
      NSITES  = 0
      NSTOT   = 0
      IATOM=0
      IMOL=0
      FNSTR	= 0.
      IST=0
      TOTCH=0.
      if(IPRINT.ge.9)then
	  write(*,*)
	  write(*,*)' Check references:'
	  write(*,*)' Iatom    Imol    Isit   Ityp     Name      Nb  '
	end if
	DO ITYP = 1,NTYPES
          NOP     = NOP    + NSPEC(ITYP)
          if(NOP.gt.NPART)stop '!!! Increase NPART in mdee.h'
          NSITES  = NSITES + NSITS(ITYP)
          IF(NSITES.GT.NS   ) STOP '!!! INCREASE NS    !'
          if(NSITS(ITYP).eq.2)FNSTR=FNSTR+2.d0*NSPEC(ITYP)
          if(NSITS(ITYP).ge.3)FNSTR=FNSTR+3.d0*NSPEC(ITYP)
          NSTOT   = NSTOT  + NSPEC(ITYP)*NSITS(ITYP)
          IF(NSTOT .GT.NTOT ) STOP '!!! INCREASE NTOT  !'
          FNSS3   = 3.D0*DFLOAT(NSITES)
          IF(NSPEC(ITYP).LE.0) STOP '!!! CHECK NSPEC !!!'
                 TFACT(ITYP) = 1.D0
          IF(LSHEJK) TFACT(ITYP) = FNSS3/(FNSS3-NRB(ITYP))
          NRTSC(ITYP) = 0
          write(*,'(a6,i4,a3,a6,3x,a5,i5,4x,a6,i4)')
     +' type ',ITYP,' - ',NAME(ITYP),' Nmol',NSPEC(ITYP),
     +'Nsites',NSITS(ITYP)
	  do IS=1,NSPEC(ITYP)
	    IMOL=IMOL+1
            ITM(IMOL)=ITYP
            do I=1,NSITS(ITYP)
	      ISITE=IST+I
	      IATOM=IATOM+1
	      NNUM(IATOM)=IMOL
              ITS(ISITE)=ITYP
	      NSITE(IATOM)=ISITE
	      ITYPE(IATOM)=ITYP
	      TOTCH	= TOTCH+CHARGE(ISITE)
	      if(IPRINT.ge.9)write(*,'(4I6,2x,2a4,a6,2x,I5)')
     +      IATOM,IMOL,ISITE,ITYP,NM(ISITE),' on ',NAME(ITYP),NRB(ITYP)
	    end do
	  end do
	  IST=IST+NSITS(ITYP)
      END DO! OF ITYP
      write(*,*)
      write(*,'(a,f8.3)')'*** Total charge            Q =  ',TOTCH
      IF(NOP.EQ.1) RHO = 0.D0
*
      PRINT "(/'*** NR OF SITES/MOL (NSITES)  ',I6)" ,NSITES
      PRINT "( '*** NR OF SITES/TOT (NSTOT)   ',I6)" ,NSTOT
      PRINT "( '*** NR OF PARTICLES (NOP)     ',I6)",NOP
	PRINT "( '*** NR OF BONDS (NRBT)        ',I6)",NRBT
*
*
*  FNST - number of degrees of freedom 
*  Obs!  total momenta is fixed! 
      FNST        = 3.d0*DFLOAT(NSTOT)-3.d0
      if(LSHEJK)FNST=FNST-NRBT
	FTIMES      = 1.D0/DFLOAT(NFREQ) 
      FNOP        = DFLOAT(NOP)
      FNSTI	= FNST-3.d0*FNOP-FNSTR+3.d0
      write(*,'(a,f6.0)')'*** NR of DEGREES of FREEDOM ',FNST
      write(*,'(a,f6.0)')'***              rotational: ',FNSTR
      write(*,'(a,f6.0)')'***              internal:   ',FNSTI
*     
	write(*,*)
      CALL IUNITS
      CALL ELRCLJ
	write(*,*)
*
*  concentrations in M/l
*     
	write(*,'(a)')'   Concentrations of different molecules:'
	do ITYP=1,NTYPES
	  CONC=NSPEC(ITYP)/(AVSNO*1.d-27*VOL)
	  write (*,'(4x,a4,f9.4,a4)')
     +  NAME(ITYP)(1:4),CONC,' M/l'  
	end do
*  Zeros
	do I=1,NRH
	  do J=1,LHIST
	    HIST(I,J)=0.
	  end do
	end do
*
*	Parameters for NVT-ensemble:
*
      IF(LNVT) THEN
        DQT	= (UNITT*1.d15/QT)**2
	  write(*,*)
        write(*,'(a,f10.3)')'*** Constant temperature MD: T = ',TRTEMP
        PRINT *,'   NOSE-NVT RELAXATION TIME      ',QT,' fs'
        LSCLT      = .FALSE.
      END IF!(LNVT...
*
*	Parameters for NPT-ensemble:
*
      IF(LNPT) THEN
        if(.not.LNVT)write(*,*)'!!! Attention: NPE-ensemble???'
        DQT	= (UNITT*1.d15/QT)**2